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Abstract 

The  overall  objective  of  the  research  was  to  advance  and  extend  methodologies  for  the  modeling 
and  simulation  of  turbulent  combustion.  Probability  density  function  (PDF)  calculations  were 
performed  of  piloted  jet  turbulent  non-premixed  flames  to  demonstrate  the  ability  of  the 
approach  to  represent  local  extinction  and  rc-ignition  and  to  characterize  the  performance  of 
different  sub-models.  PDF  calculations  were  performed  of  lifted  flames  in  vitiated  co-flows,  and 
the  results  were  analyzed  to  show  that  the  stabilization  mechanism  is  auto-ignition.  Time¬ 
averaging  strategies  were  shown  to  be  very  effective  in  reducing  bias  in  the  Monte  Carlo 
methods  used  to  solve  the  PDF  equations.  Strategies  were  developed  for  implementing 
combustion  chemistry  on  parallel  computers,  resulting  in  speed-ups  of  up  to  a  factor  of  five 
(relative  to  the  straightforward  purcly-local  strategy).  The  method  was  shown  to  scale  up  to 
4096  cores.  Algorithms  and  codes  were  developed  to  implement  the  combined  methodology  of 
large-eddy  simulation  (LES)  and  filtered  density  function  (FDF).  Second-order  schemes  were 
developed  and  demonstrated  for  particle  tracking,  which  (in  contrast  to  standard  methods)  also 
accurately  represent  the  evolution  of  the  particle  number  density,  and  second-order  accurate 
splitting  schemes  were  developed  for  the  stochastic  differential  equations  that  arise  in  LES/FDF. 


20090429219 


1 


AFRL-SR-AR-TR-09-0 1 32 


REPORT  DOCUMENTATION  PAGE 

Public  reporting  burden  for  this  coiiection  of  information  is  estimated  to  average  1  hour  per  response,  including  tha  tima  for  raviawing  a 

data  needed,  and  completing  and  reviewing  this  collection  erf  information.  Send  comments  regarding  this  burden  estimata  or  any  other  aspeci  ui  um3  -  g 

this  burden  to  Dapartmant  of  Defanse  Washington  Headquartarj  Services.  Directorata  for  information  Operations  and  Raports  (0704-0188).  1215  Jaffarson  Davis  Highway.  Suite  1204,  Arlington.  VA  2^u2- 
4302.  Respondtnts  should  be  awara  that  notwithstanding  any  othar  provision  of  law  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

23-2-2009  Final  Performance  Report 

3.  DATES  COVERED  (From  -  To) 
1-1-2006-31-11-2008 

4.  TITLE  AND  SUBTITLE 

(U)  PDF  Modelling  of  Turbulent  Combustion 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

FA9550-06-, *1-0048 

5c.  PROGRAM  ELEMENT  NUMBER 

61102F 

6.  AUTHOR(S) 

Stephen  B.  Pope 

5d.  PROJECT  NUMBER 

2308 

5e.  TASK  NUMBER 

BX 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Cornell  University 

Upson  Hall 

Ithaca  NY  14853 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR/NA 

875  Randolph  St 

Suite  325,  Room  3112 

Arlington  VA  22203-1768 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

The  overall  objective  of  the  research  was  to  advance  and  extend  methodologies  for  the  modeling  and  simulation  of  turbulent  combustion. 
Probability  density  function  (PDF)  calculations  were  performed  of  piloted  jet  turbulent  non-premixed  flames  to  demonstrate  the  ability  of 
the  approach  to  represent  local  extinction  and  re-ignition  and  to  characterize  the  performance  of  different  sub-models.  PDF  calculations 
were  performed  of  lifted  flames  in  vitiated  co-flows,  and  the  results  were  analyzed  to  show  that  the  stabilization  mechanism  is  auto¬ 
ignition.  Time-averaging  strategies  were  shown  to  be  very  effective  in  reducing  bias  in  the  Monte  Carlo  methods  used  to  solve  the  PDF 
equations.  Strategies  were  developed  for  implementing  combustion  chemistry  on  parallel  computers,  resulting  in  speed-ups  of  up  to  a 
factor  of  five  (relative  to  the  straightforward  purely-local  strategy).  The  method  was  shown  to  scale  up  to  4096  cores.  Algorithms  and 
codes  were  developed  to  implement  the  combined  methodology  of  large-eddy  simulation  (LES)  and  filtered  density  function  (FDF). 
Second-order  schemes  were  developed  and  demonstrated  for  particle  tracking,  which  (in  contrast  to  standard  methods)  also  accurately 
represent  the  evolution  of  the  particle  number  density,  and  second-order  accurate  splitting  schemes  were  developed  for  the  stochastic 
differential  equations  that  arise  in  LES/FDF. 

15.  SUBJECT  TERMS 

Propulsion,  combustion,  PDF  methods 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

UL 

18.  NUMBER 
OF  PAGES 

62 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Julian  M.  Tishkoff 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

19b.  TELEPHONE  NUMBER  (include area 
code) 

(703)  696  8478 

Standard  Form  298  {Rev.  8-98) 

Prescribed  by  ANSI  Std  239.18 


CONTENT 


EXECUTIVE  SUMMARY . 3 

RANS/PDF  CALCULATIONS  OF  TURBULENT  FLAMES . 4 

PDF  CALCULATIONS  OF  PILOTED  JET  FLAMES . 4 

PDF  CALCULATIONS  OF  LIFTED  TURBULENT  FLAMES . 7 

TIME-AVERAGING  STRATEGIES  IN  PDF  METHODS . 8 

EFFICIENT  PARALLEL  IMPLEMENTATION  OF  IN  SITU  ADAPTIVE  TABULATION . 9 

ALGORITHMS  FOR  LES/FDF . 10 

ACCURATE  PARTICLE  TRACKING . 10 

TIME  SPLITTING  SCHEMES . 14 

PUBLICATIONS . 16 

OTHER  REFERENCES . 16 

PERSONNEL  SUPPORTED . 18 

DEGREES  GRANTED . 18 

MEETING  PARTICIPATION . 18 

INTERACTION  WITH  AFRL . 19 

TECHNOLOGY  TRANSITIONS  AND  TRANSFERS . 20 

HONORS  AND  AWARDS . 20 

INVENTIONS  AND  PATENTS . 20 


2 


EXECUTIVE  SUMMARY 


In  both  space  and  aircraft  applications,  the  design  of  combustors  in  propulsion  systems  remains  a 
significant  technical  challenge.  Given  the  cost,  difficulty,  and  time  consumed  in  experimental 
testing,  it  is  recognized  well  that  computer  modeling  is  essential  to  exploring  different  design 
concepts  and  to  reducing  the  cost  and  time  of  the  design  cycle.  While  many  phenomena  may  be 
involved  -  sprays,  radiation,  combustion  dynamics,  etc.  -  a  central  problem  is  that  of  modeling 
turbulent-chemistry  interactions  in  turbulent  combustion.  The  PDF  approach  to  turbulent 
combustion  has  the  advantages  of  fully  representing  the  turbulent  fluctuations  of  species  and 
temperature  and  of  allowing  realistic  combustion  chemistry  to  be  implemented  (e.g.,  of  order  50 
species).  This  methodology  also  is  being  applied  in  conjunction  with  large-eddy  simulations,  in 
which  case  it  is  referred  to  as  LES/FDF. 

The  research  performed  and  findings  obtained  are  summarized  here  and  given  in  more  detail  in 
the  subsequent  sections  and  in  the  listed  publications. 

Turbulent  flames  were  studied  using  probability  density  function  (PDF)  methods.  The  influence 
of  the  mixing  sub-models  was  investigated  for  a  series  of  piloted  jet  non-premixed  flames.  Only 
one  of  the  frequently  used  models  (EMST,  Subramaniam  &  Pope  1998)  was  capable  of 
predicting  the  correct  levels  of  both  the  degree  of  local  extinction  and  the  level  of  turbulent 
fluctuations. 

A  Lagrangian  study  of  the  performance  of  different  turbulent  mixing  models  was  performed  for 
lifted  jet  flames  in  vitiated  co-flows.  This  study  revealed  that  the  stabilization  mechanism  in 
these  flames  is  autoignition,  following  the  turbulent  mixing  of  the  jet  and  coflow.  In  contrast  to 
conventional  lifted  flames,  no  diffusive  transport  against  the  flow  was  needed  for  stabilization. 

Computational  algorithms  were  developed  to  accelerate  parallel  combustion  chemistry 
calculations.  For  serial  calculations,  the  method  of  in  situ  adaptive  tabulation  (ISAT)  was  able  to 
accelerate  chemistry  calculations  by  three  orders  of  magnitude.  For  parallel  computations  there 
were  serious  load-balancing  problems.  Several  different  algorithms  were  implemented  and 
evaluated,  which  partially  resolved  this  problem.  Speed-ups  of  up  to  five  were  obtained,  and 
scalability  to  at  least  4096  cores. 

Advances  have  been  made  in  the  development  of  accurate  and  efficient  computational  algorithms 
for  PDF  and  LES/FDF  methods.  Time-averaging  strategies  in  PDF  methods  were  shown  to  be 
very  effective  in  reducing  the  bias  error.  New  interpolation  schemes  and  a  modified  Runge 
Kutta  scheme  were  shown  to  yield  superior  performance  for  particle  tracking  in  LES/FDF. 
These  schemes  are  second-order  accurate  in  time  and  space  and  cause  the  particle  number 
density  to  evolve  accurately.  Second-order  accurate  splitting  schemes  were  developed  for  the 
stochastic  differential  equations  that  arise  in  LES/FDF. 

The  grant  provided  support  for  the  PI  and  for  three  Ph.D.  students:  Liuyan  Lu,  Haifcng  Wang 
and  Pavel  Popov.  Lu  received  her  Ph.D.  in  2008. 
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RANS/PDF  CALCULATIONS  OF  TURBULENT  FLAMES 


PDF  CALCULATIONS  OF  PILOTED  JET  FLAMES 

The  Barlow  &  Frank  (1998)  piloted  jet  flames  are  recognized  well  as  providing  an  excellent  test 
of  turbulent  combustion  models,  in  particular  of  their  ability  to  describe  local  extinction  and  re¬ 
ignition.  There  is  a  sequence  of  flames  A-F,  with  D,  E,  and  F  being  fully  turbulent,  with 
increasing  amounts  of  local  extinction.  While  there  are  many  successful  modeling  studies  of 
flame  D  —  which  has  little  local  extinction  -  there  are  far  fewer  modeling  studies  for  the  more 
challenging  flames  E  and  F. 

PDF  methods  previously  have  been  applied  successfully  to  these  flames  by  the  Cornell  group 
(Xu  &  Pope,  2000,  Tang  et  al.,  2000)  and  the  Imperial  College  group  (Lindstedt  et  ah,  2000). 
Both  groups  reported  accurate  calculations  of  these  flames  but  with  different  sub-models.  The 
former  calculations  use  an  augmented  reduced  mechanism  based  on  the  GRI  mechanism  and  the 
EMST  mixing  model  with  model  constant  C $  =  1.5.  The  latter  calculations  use  Lindstedt’s 
reduced  mechanism  and  the  modified  Curl  mixing  model  with  model  constant  Q>=  2.3.  Since 
the  two  groups  use  different  mixing  models  and  different  chemical  mechanisms,  questions  arise 
as  to  the  dependence  and  sensitivity  of  the  calculations  to  these  ingredients.  Some  progress  has 
been  made  in  understanding  the  relative  behavior  of  the  different  mixing  models  (see,  e.g.,  Ren 
&  Pope  2004).  We  now  have  concluded  a  series  of  PDF  calculations  addressing  these  issues. 
Cao  &  Pope  (2005)  describe  joint  PDF  calculations  of  these  flames  using  7  different  chemical 
mechanisms  for  methane,  ranging  from  a  5-step  reduced  mechanism  to  the  53-species  GRI  3.0 
mechanism.  The  results  show  that  the  GRI  mechanisms  (and  the  augmented  reduced 
mechanisms  based  on  them)  are  capable  of  accurately  representing  the  observed  local  extinction 
and  rc-ignition.  In  contrast,  Ci  skeletal  mechanisms  and  5-step  reduced  mechanisms  prove  to  be 
inaccurate.  (Lindstedt’s  mechanism  was  not  included  in  the  study  due  to  its  lack  of  availability 
at  the  time.) 

This  work  has  been  followed  up  by  a  study  (Cao  et  al.  2006)  of  the  influence  of  the  turbulent 
mixing  models  used  in  the  PDF  model  calculations.  The  principal  finding  is  that  (with  a 
different  value  of  the  mixing-model  constant)  all  three  commonly  used  models  (IEM,  MC  and 
EMST)  are  capable  of  yielding  the  observed  level  of  local  extinction,  but  only  the  EMST  model 
produces,  simultaneously,  the  observed  level  of  mixture  fraction  fluctuations. 

After  Lindstedt’s  mechanisms  became  available  publicly  on  the  TNF  workhop  website 
(http://public.ca.sandia.gov/TNF/chemistry.html),  we  further  investigated  the  apparent 
discrepancies  in  previous  PDF  model  calculations  of  the  Barlow  &  Frank  flames  (Wang  and 
Pope,  2007).  A  comparison  is  made  between  non-premixed  flame  calculations  using  the  GRI  3.0 
and  Lindstedt  detailed  mechanisms  and  with  different  mixing  models.  Additional  insights  arc 
provided  by  autoignition  and  non-premixed  strained  laminar  flame  calculations.  Figure  1  shows 
the  ignition  delay  times  given  by  different  mechanisms  for  a  stoichiometric  methane/air  mixture 
as  a  function  of  the  initial  temperature.  For  temperatures  above  1,400K,  the  Lindstedt 
mechanism  exhibits  a  shorter  ignition  delay  time.  A  series  of  calculations  was  performed  of 
strained  non-premixed  laminar  flames  with  the  same  fuel  and  oxidant  used  in  the  Barlow  & 
Frank  flames  (i.e.,  the  fuel-air  ratio  is  1:3  methanc/air,  and  the  oxidant  is  air).  Figure  2  shows 
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the  peaks  of  various  quantities  in  these  flames  as  functions  of  the  imposed  strain  rate.  While  the 
Lindstedt  and  GRI  mechanisms  yield  comparable  properties  at  low  strain  rates,  the  Lindstcdt 
mechanism  exhibits  a  greater  resilience  to  extinction. 
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Figure  1  Ignition  delay  times  (IDTs)  of  different  mechanisms  at  different  initial  temperatures  T. 
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Figure  2  The  peak  values  of  the  temperature  and  the  mass  fractions  of  CO2,  H2O,  CO,  H2,  and 
OH  against  the  nominal  strain  rate  a  in  the  opposed  laminar  jet  non-premixed  flames  computed 
with  different  reaction  mechanisms. 
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A  set  of  28  PDF  calculations  of  the  Barlow  &  Frank  flames  have  been  performed,  with  both 
chemical  mechanisms,  three  mixing  models,  and  a  range  of  values  of  the  mixing  model  constant 
Cfr  Figure  3  shows  some  of  the  results  on  “burning  indices”  (Bis)  at  an  axial  location  of  7.5  ftiel 
jet  diameters.  (The  burning  index  is  zero  for  complete  extinction  and  one  for  complete 
combustion.)  There  are  two  fundamental  observations  from  these  results.  First,  for  a  given 
mixing  model  and  value  of  the  Lindstedt  mechanism  uniformly  yields  higher  values  of  BI 
than  does  the  GRI  mechanism.  This  result  is  consistent  with  the  previously  observed  shorter 
ignition  delay  time  and  higher  extinction  strain  rate  of  the  Lindstedt  mechanism.  Second,  the 
two  calculations  closest  to  those  of  Xu  &  Pope  (2000)  and  Lindstedt  et  al.  (2000)  (i.e.,  GRI, 
EMST,  =  1.5  and  Lindstedt,  modified  Curl,  C,p  =  2.3)  yield  very  similar  values  of  BI,  which 
are  generally  in  agreement  with  the  experimental  data  (though  less  so  for  H2O  and  OH).  This 
result,  then,  provides  a  full  explanation  for  the  similar  calculations  in  2000  using  different 
submodels. 


Figure  3  The  burning  indices  of  the  temperature  and  the  mass  fractions  of  CO2,  H2O,  CO,  H2, 
and  OH  against  C^>  at  the  location  of  x/D  =  7.5  in  the  Sandia  flame.  The  two  letters  in  the  legend 
identify  the  combination  of  the  mixing  models  (E:  ESMT,  I:  IEM,  and  M:  modified  Curl)  and  the 
reaction  mechanisms  (G:  GRI-Mech  3.0,  L:  Lindstedt  mechanism). 

The  present  work  reveals  significant  differences  in  the  predictions  of  the  two  mechanisms  but 
does  not  determine  which  (if  either)  is  accurate  under  the  conditions  investigated.  There  are 
experimental  data  for  non-premixed  laminar  flame  with  the  flame  E  fuel  (mcthane/air  mixture 
Barlow  et  al.  2001),  but  these  are  at  very  low  strain  rates.  It  would  be  extremely  valuable  to 
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investigate  experimentally  the  properties  of  strained  laminar  flames  of  this  fuel  close  to  the 
extinction  strain  rate. 

It  was  claimed  later  by  Lindstedt  (private  communication)  that  the  Lindstcdt  mechanism 
provided  on  the  TNF  web  site  and  used  in  the  above  investigation  includes  incorrect 
specifications  of  reaction  parameters.  As  a  consequence,  the  above  conclusions  on  the 
comparison  between  the  GRI  and  the  Lindstcdt  mechanisms  arc  compromised.  The  difference 
between  the  two  sets  of  PDF  calculations  in  2000  unfortunately  remains  unexplained,  and  further 
investigation  is  needed. 

PDF  CALCULATIONS  OF  LIFTED  TURBULENT  FLAMES 

In  laboratory  experiments  on  turbulent  flames,  for  obvious  reasons,  air  at  atmospheric  conditions 
generally  is  used  as  the  oxidant.  In  practical  applications,  however,  the  re-circulating  flows  used 
for  flame  stabilization  generally  lead  to  some  mixing  between  the  air  stream  and  hot  combustion 
products,  and  in  thrust  augmentors  the  oxidant  stream  is  simply  the  lean  combustion  products 
from  the  turbine. 
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Figure  4  Axial  profiles  (off  the  axis  at  r=15mm)  of  mean  mixture  fraction,  temperature,  and 
mass  fractions  of  CH4,  O2,  CO2,  H2O,  CO,  H2  and  OH.  Symbols:  Measurement  by  Cabra  et  al. 
2005;  Solid  lines:  Joint  PDF  calculations  with  ISAT  error  tolerance  fetol  =5x10  4 ;  Dashed  lines: 

^etoi  =1x10'"’;  Dash-dotted  lines:  fetol  =  lxl0~5;  Dotted  lines:  fetol  =  5xl0~6. 


Motivated  by  these  observations,  a  series  of  experiments  has  been  performed  (at  Berkeley, 
Sandia,  and  Sydney)  on  various  jet  flames  in  vitiated  co-flows.  Previously  we  have  performed 
two  PDF  studies  of  the  lifted  hydrogen  flames  studied  experimentally  by  Cabra  et  al.  (2002). 
These  two  studies  are  based  on  the  composition  PDF  method  incorporated  in  Fluent  (Masri  et  al. 
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2004)  and  on  the  velocity-frequency-composition  joint  PDF  method  (Cao,  Pope  &  Masri  2005). 
In  both  cases,  detailed  9-specics  mechanisms  for  hydrogen  arc  used. 

In  this  study,  we  investigate  the  lifted  methane  flames  studied  by  Cabra  et  al.  (2002,  2005).  The 
calculations  are  based  on  the  joint  PDF  method  used  previously,  with  the  modified  Curl  mixing 
model  and  the  ARM1  mechanism.  Figure  4  shows  the  comparison  between  experimental  data 
and  PDF  results  with  different  values  of  the  ISAT  (Pope,  1997)  error  toleranee.  In  order  to 
obtain  accurate  PDF  results  of  the  lifted  methane  flames,  the  ISAT  error  toleranee  needs  to  be 
decreased  to  at  least  £etol  =  lxlCT5  for  these  flames.  Similar  to  the  lifted  hydrogen  flames,  strong 
sensitivity  of  the  results  to  the  coflow  temperature  is  observed  in  the  calculations. 

TIME-AVERAGING  STRATEGIES  IN  PDF  METHODS 

Computationally,  PDF  methods  usually  are  implemented  using  Lagrangian  particle  Monte  Carlo 
methods.  With  N  being  the  number  of  particles  used,  these  methods  incur  a  statistical  error  of 
order  N  VJ  and  a  bias  error  of  order  N  1 .  At  first  sight,  the  statistical  error  appears  dominant; 
however,  unlike  the  bias,  the  statistical  error  has  zero  mean,  and  for  statistieally-stationary 
problems  it  can  be  reduced  at  will  by  time  averaging.  Hence,  in  PDF  methods,  bias  can  be  the 
more  serious  source  of  numerical  error. 
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Figure  5  The  means  of  the  axial  velocity  U,  the  turbulence  frequency  a> ,  the  Reynolds  stresses 
uu  and  uv,  and  the  means  and  rms  of  the  mixture  fraction  and  temperature  against  N~'  at  the 

location  of  (15D,  ID)  in  the  Cabra  lifted  H2/N2  jet  flame.  Circles,  with  time-averaging; 
diamonds,  without  time-averaging 


Wang  &  Pope  (2008b)  show  how  the  bias  can  be  reduced  dramatically  in  the  veloeity- 
composition  PDF  method.  In  this  approach  the  PDF  considered  is  for  the  fluctuating  component 
of  velocity  (and  composition).  A  correction  is  applied  to  enforce  the  consistency  condition  that 
the  mean  of  the  velocity  fluctuation  is  zero.  This  correction  can  be  applied  based  either  on  the 
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instantaneous  mean  of  the  fluctuating  particle  velocities  or  on  the  time  average  of  this  quantity. 
As  demonstrated,  the  use  of  time  averaging  all  but  eliminates  the  substantial  bias  which  is 
otherwise  incurred. 

EFFICIENT  PARALLEL  IMPLEMENTATION  OF  IN  SITU  ADAPTIVE 
TABULATION 

The  PDF  method  calculations  presented  above  demonstrate  current  capabilities,  e.g.,  using  a 
mechanism  with  order  50  species.  Such  calculations  are  possible  because  of  the  ISAT  algorithm 
(Pope  1997)  and  the  availability  of  parallel  computers.  As  we  work  towards  combining  PDF 
methods  with  large  eddy  simulation  (LES),  larger  parallel  clusters  are  needed,  and  the  efficient 
parallel  implementation  of  ISAT  becomes  crucial. 

If  there  are  P  processors,  then  the  load  associated  with  solving  the  LES  equations  is  balanced 
well  if  the  solution  domain  is  partitioned  into  P  sub-domains,  each  containing  approximately  the 
same  number  of  cells.  Further,  with  such  a  decomposition,  the  computational  particles  used  in 
the  PDF  algorithm  are  distributed  approximately  equally  among  the  processors;  hence,  there  is 
good  load  balancing  of  particle  tracking  and  of  mixing.  The  load  associated  with  reaction  (i.c., 
incrementing  the  particle  compositions  due  to  reaction  over  the  time  step)  that  can  dominate  the 
CPU  time  can  be  poorly  balanced.  For  example,  the  particles  on  one  processor  may  be  inert 
(e.g.,  cold  air)  and  hence  require  negligible  time  to  evaluate  their  reaction,  whereas  the  particles 
on  another  processor  may  be  highly  reactive.  Even  using  ISAT,  reactive  particles  take  much 
longer  to  treat,  especially  if  retrieving  from  a  full  table  is  not  possible;  hence,  a  direct  integration 
of  the  stiff  ODEs  is  required. 


Figure  6  Contour  plot  of  mean  temperature  in  PDF  calculations  of  a  lifted  hydrogen  jet  flame, 
showing  the  sub-domains  used  in  the  Fluent  parallel  computation. 


Following  the  preliminary  work  of  Lu  et  al.  (2005),  we  are  continuing  to  develop  and  evaluate 
parallel  strategies  for  implementing  ISAT.  One  of  the  test  problems  considered  is  the 
composition  PDF  method  in  Fluent  applied  to  the  lifted  hydrogen  jet  flame  described  above. 
Figure  6  shows  the  solution  domain  and  its  decomposition  into  the  eight  sub-domains  used  in  the 
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parallel  computation  (using  8  processors).  Figure  7  shows  the  CPU  and  wall-clock  times  per 
particle  step  required  by  three  different  parallel  strategies.  The  first  is  PLP  (purely  local 
processing)  in  which  no  message  passing  is  performed,  but  each  processor  invokes  ISAT 
independently  for  all  of  the  particles  in  its  sub-domain.  As  may  be  seen  from  the  CPU  times, 
there  is  considerable  load  imbalance  with  processor  6  requiring  significantly  more  time.  The 
corresponding  sub-domain  is  the  base  of  the  flame,  which  involves  the  most  challenging  ignition 
chemistry.  The  two  (non-trivial)  parallel  strategies  lead  to  speed-ups  (relative  to  PLP)  of  about 
2.5  and  5. 
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Figure  7  Average  wall  clock  and  CPU  times  (per  particle  per  reaction  fractional  step)  for  parallel 
Fluent  computations  of  the  lifted  hydrogen  flame,  with  different  parallel  ISAT  strategics. 

An  adaptive  pairing  algorithm  (Lu  et  al.,  2008)  balances  the  cost  of  direct  evaluations  and  the 
communication  time  for  distributing  particles  among  processors.  This  algorithm  performs 
groupings  of  processors  based  on  the  probability  of  retrieving  from  those  ISAT  tables,  the  cost  of 
distributing  the  particles,  and  the  cost  of  performing  direct  evaluations.  Particles  are  distributed 
among  the  processors  in  the  groups  using  the  aforementioned  parallel  strategies.  Partially-stirred 
reactor  test  cases  have  shown  functionality  and  scalability  on  up  to  4,096  cores  on  Ranger  at  the 
Texas  Advanced  Computing  Center. 

ALGORITHMS  FOR  LES/FDF 

The  bulk  of  the  research  effort  has  been  devoted  to  developing  algorithms  and  code  for  the 
application  of  the  LES/FDF  approach  to  turbulent  flames.  Here  we  report  two  important 
algorithms  in  FDF,  namely  accurate  tracking  of  particles  through  a  given  velocity/diffusivity 
field  and  weak  second-order  time  splitting  schemes  for  solving  the  particle  stochastic  differential 
equations.  The  development  and  testing  of  the  in-house  FDF  code  is  close  to  being  finished  and 
will  be  brought  to  the  full  LES/FDF  studies  of  turbulent  combustion  soon. 

ACCURATE  PARTICLE  TRACKING 
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One  important  issue  in  LES/FDF  numerical  methods  is  the  consistency  between  particle  number 
density,  which  can  be  thought  of  as  the  expected  number  of  particles  in  a  unit  volume,  and  the 
LES  filtered  density.  To  ensure  the  validity  of  the  LES/FDF  methodology,  these  two  redundant 
forms  of  density  must  be  equal.  It  is  because  of  this  fact  that  particle  tracking  schemes  that 
preserve  the  consistency  between  these  two  fields,  either  explicitly  or  implicitly,  have  received 
much  attention  in  recent  years. 


To  illustrate  the  issues  and  advances  made,  wc  consider  the  simplest  case  of  fluid  particles  that 
move  according  to  an  incompressible  velocity  field,  U(x,t),  with  no  turbulent  diffusivity.  This 

problem,  albeit  being  a  simplification  from  the  general  stochastic  differential  equation  that 
governs  particle  transport  (see  Eqs.  (6)  and  (7)  below),  is  important  in  the  LES/FDF  context, 
where  a  significant  part  of  the  turbulent  transport  is  carried  out  by  the  resolved  velocity  field. 
The  particle  position  X(t)  evolves  by  the  ODE 


^P=u(x(,),,) 


(1) 


Let  q(x,t)  denote  the  particle  number  density,  and  let  this  be  initially  uniform,  ^(x,0)  -q0 .  Just 
like  the  fluid  density,  the  particle  number  density  satisfies  the  conservation  equation. 

(2) 


§+''’•(«*')  =  0. 


which  can  be  re-expressed  as 


\ 


\nq  =  -VU 


-  +  U-V 

ydt 

Thus,  given  V  •  U  =  0 ,  it  follows  that  an  initially  uniform  distribution  remains  uniform. 


(3) 


The  above  results  pertain  to  an  exact  implementation.  In  practice  the  velocity  field  is  represented 
discretely  on  a  grid,  and  Eq.  (1)  is  integrated  numerically,  with  U(  X(t),t)  being  approximated 

by  some  interpolation.  Consider  2D  for  simplicity,  and  suppose  that  (as  is  usually  the  case)  bi¬ 
linear  interpolation  is  used.  This  interpolation  scheme  implies  that  the  velocity  experienced  by 
the  particles  inside  each  grid  cell  is  of  the  form 

u  =  u0+uxx+uyy  +  uxyxy, 

v  =  vo  +  V*x  +  Vvy  +  Vxyxy. 

Evidently,  the  implied  velocity  divergence  is 

_  du  dv 

u=!k+ty=Ux+Uiyy+Vy+ VxyX'  (5) 

Which,  in  general  (i.e.,  urt  *  0,vxi  *  0),  cannot  be  zero. 


(4) 


Figure  8  shows  the  result  of  a  simple  test  in  which  particles  (initially  uniform)  are  tracked 
through  an  incompressible  velocity  field  using  bilinear  interpolation.  Clearly,  the  non-zero 
velocity  divergence  experienced  by  the  particles  leads  to  a  decidedly  non-uniform  distribution. 

To  overcome  this  problem,  McDermott  &  Pope  (2008)  devised  a  “parabolic  edge  reconstruction 
method”  of  interpolation  (PERM)  that  (for  the  case  considered)  yields  an  interpolated  velocity 
field  of  zero  divergence;  thus,  if  the  fluid  particles  are  advanced  accurately  in  time,  the  uniform 
initial  distribution  is  preserved. 


0  1  2  3  4  5  6 

Figure  8  Left:  Exact  solution  for  an  incompressible  2D  flow.  Right:  numerical  solution  for  the 
same  flow,  using  bilinear  velocity  interpolation 


RK2,  one  flow-through  MRK2,  one  flow-through 

Figure  9  Particles  tracked  through  a  2D  incompressible  flow  field  using  PERM  interpolation  and 
(a)  (left)  second-order  Runge-Kutta  (b)  (right)  MRK2.  The  Courant  number  is  1 .0.  In  both  cases 
the  non-uniformity  decreases  as  the  Courant  number  decreases. 

The  PERM  velocity  field,  by  construction,  has  excellent  divergence  properties;  however,  it  also 
results  in  small  velocity  discontinuities  between  cells  (that  decrease  quadratically  as  the  cell  size 
is  decreased).  It  follows  that  a  standard  numerical  scheme  for  integrating  Eq.  (1)  reverts  to  first- 
order  accuracy  in  the  time  step  At  (because  of  the  discontinuities).  To  remedy  this  problem, 
Popov  et  al.  (2008)  have  devised  a  modified  Runge-Kutta  scheme  (MRK2)  that  restores  second- 
order  accuracy  even  with  the  type  of  discontinuity  caused  by  PERM.  Figure  9  shows  the  particle 
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distribution  for  a  test  case  using  (a)  standard  second-order  Runge-Kutta,  and  (b)  MRK2.  As  may 
be  seen,  MRK2  results  in  a  much  more  uniform  number  density. 


Figure  10  Left:  Solution  for  an  incompressible  flow  on  a  circular  domain,  using  bilinear  velocity 
interpolation  in  polar  coordinates.  Right:  Solution  for  the  same  flow  using  the  PPERM  velocity 
reconstruction 


Figure  1 1  Overall  numerical  error  of  the  particle  tracking  scheme,  composed  of  the  velocity  and 
scalar  interpolation  schemes,  and  the  SDE  time  integration  schemes  developed  by  Cao  and  Pope 
(CP,  red  dots),  and  Tocino  and  Vigo-Aguiar  (TVA,  blue  dots).  The  vertical  bars  represent  95% 
confidence  intervals. 

Because  a  cylindrical  grid  is  the  most  natural  choice  for  many  turbulent  combustion  simulations, 
recent  efforts  have  included  the  adaption  of  the  PERM  velocity  reconstruction  algorithm  (that  is 
designed  for  Cartesian  grids)  to  cylindrical  grids.  Because  such  a  grid  is  inherently  singular  at  its 
axis,  care  must  be  taken  to  avoid  singularities  in  the  reconstructed  velocity,  or  else  a  decrease  in 
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the  order  of  accuracy  of  the  scheme  is  suffered,  close  to  the  grid  axis.  The  new,  polar  version  of 
PERM  (called  PPERM),  avoids  these  complications  and  yields  an  interpolated  velocity  field 
whose  accuracy  and  regularity  properties  are  identical  to  those  of  PERM.  The  performance  of 
the  PPERM  scheme  is  depicted  on  Figure  10,  which  shows  results  from  a  test  analogous  to  the 
one  illustrated  on  Figure  9,  performed  on  a  circular  domain. 


Additional  efforts  in  particle  tracking  have  included  testing  of  the  overall  accuracy  of  the 
combination  of  interpolation  and  time  advancement  schemes  for  the  general,  stochastic  particle 
tracking  problem.  In  the  majority  of  composition  LES/FDF  methods,  the  particle  positions 
evolve  by  the  stochastic  differential  equation 


<*(')  = 


vr 


(6) 


u+— 

V  p  ) 

where  T  is  the  effective  diffusivity,  p  is  the  LES  filtered  density,  the  final  term  in  Eq.  (6)  has 
the  standard  meaning  of  integration  with  respect  to  a  Wiener  process.  Velocity  is  interpolated 
with  the  PPERM  reconstruction  scheme,  and  the  diffusivity  and  density  fields  are  interpolated 
with  a  third-order  accurate  scheme,  called  the  multilinear  gradients  (MLG)  scheme.  Two 
alternative  SDE  time  advancement  schemes,  developed  respectively  by  Cao  and  Pope  (2003)  and 
Tocino  and  Vigo- Aguiar  (2002),  have  been  tested,  and  results  have  been  reported  by  Popov  and 
Pope  (2008).  As  can  be  seen  on  Figure  1 1,  the  overall  scheme  attains  the  desired  second-order 
accuracy  using  either  of  these  SDE  schemes. 


TIME  SPLITTING  SCHEMES 


The  Lagrangian  Monte  Carlo  particle  methods  (Pope,  1985)  are  the  most  efficient  numerical 
algorithm  to  solve  PDF/FDF  transport  equations.  In  these  methods,  the  continuous  PDF  is 
discretized  by  a  finite  number  of  nominal  particles,  and  each  particle  is  governed  by  a  system  of 
stochastic  differential  equations  (SDE)  as  follows  (including  an  Ito  SDE  for  particle  position  and 
a  scalar  random  equation)  to  describe  the  underlying  physical  and  chemical  processes. 

fdX(/)  =  D(X(/),/)df  +  6(X(/),/)dW(/), 
l</#(/)  =  A(X(/),#(/),#)<*, 

where  X(/)  and  <j)(t )  are  particle  position  and  scalars  (e.g.,  species  mass  fractions  and  enthalpy), 
D(X(/),/)  and  />(X(/),/)  denote  drift  and  diffusion  in  physical  space,  and  A(X(/), 
denotes  rate  of  change  of  scalars  (e.g.,  production  rate  of  species  or  heat  release  rate).  Eq.  (6)  is 
a  specific  version  of  the  Ito  SDE  in  Eq.  (7),  with  D(X(/),/)  =  (U  + Vr/p)  and  b(X(t),t)  =  y[2 T/p . 

To  solve  the  coupled  SDE  system  (7)  accurately,  special  care  is  needed.  All  the  well  developed 
high-order  ODE  schemes  are  not  applicable  to  this  system.  No  second-order  time  splitting 
schemes  have  been  developed  previously  for  this  system  and  have  been  applied  in  the 
RANS/PDF  or  LES/FDF  practice. 

A  class  of  weak  second-order  time  splitting  schemes  is  developed  in  Wang  el  al.  (2008).  Four 
primary  contributions  are  made  in  this  work: 

I.  Establish  that  the  coefficients  in  (7)  can  be  frozen  at  the  mid-time  while  preserving 
second-order  accuracy.  The  use  of  frozen  coefficients  avoids  any  interpolation  of  the 
coefficients  in  time  when  using  different  SDE  schemes,  and  greatly  simplifies  the 
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solution  algorithms.  Moreover,  it  guarantees  the  applicability  of  many  autonomous  SDE 
schemes  to  (7). 

II.  Examine  the  performance  of  three  existing  schemes  for  integrating  the  SDE  for  particle 
position:  the  mid-point  scheme  of  Cao  and  Pope  (2003);  the  predictor-corrector  scheme 
of  Kloeden  and  Platen  (1992);  and  the  Runge-Kutta  scheme  of  Tocino  and  Vigo-Aguiar 
(2002).  The  latter  two  are  derivative-free  schemes,  and  the  Kloeden  and  Platen  scheme  is 
for  autonomous  SDEs  only. 

III.  Develop  and  evaluate  different  time  splitting  schemes  (that  treat  particle  motion,  reaction, 
and  mixing  on  different  sub-steps).  The  second-order  accuracy  of  the  schemes  is  verified 
by  the  test  cases. 

IV.  Develop  the  method  of  manufactured  solutions  (MMS)  (Roache,  2002)  to  assess  the 
convergence  of  Monte  Carlo  particle  methods.  This  method  is  a  general  way  of 
generating  test  cases  with  analytical  solutions.  The  MMS  method  is  used  to  demonstrate 
the  convergence  and  order  of  accuracy  of  the  developed  splitting  schemes. 


Figure  12  The  convergence  of  the  errors  of  scalar  mean  e-  and  scalar  square  mean  e -  against 

the  time  step  At  for  splitting  of  type  Mix-Tran-React-Tran-Mix  combined  with  the  Cao  and  Pope 
scheme  (circles),  with  the  Cao  and  Pope  scheme  with  frozen  coefficients  (diamonds),  with  the 
Tocino  and  Vigo-Aguiar  scheme  (squares),  with  the  Tocino  and  Vigo-Aguiar  scheme  with 
frozen  coefficients  (down  triangles),  and  with  the  Kloeden  and  Platen  scheme  with  frozen 
coefficients  (left  triangles).  (The  error  bars  indicate  95%  confidence  intervals.) 


Different  splitting  schemes  with  or  without  freezing  the  coefficients  are  compared.  In  general, 
freezing  the  coefficients  reduces  the  numerical  errors.  Otherwise  no  significant  differences  are 
observed  in  the  performance  of  the  different  SDE  schemes  and  splitting  schemes.  Figure  12 
shows  the  comparison  of  one  type  of  splitting  (Mix-Tran-React-Tran-Mix,  which  means  taking 
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half  step  of  mixing  first,  then  half  step  of  transport,  and  full  step  of  reaction,  then  another  half 
step  of  transport  followed  by  another  half  step  of  mixing)  with  different  SDE  schemes  and  with 
or  without  freezing  the  coefficients.  All  the  schemes  shown  on  the  plot  demonstrate  second- 
order  convergence  with  respect  to  time  step  size,  and  the  difference  among  the  different  schemes 
is  small. 

PUBLICATIONS 

The  following  papers  stemming  from  the  current  AFOSR  program  were  written  or  published  in 
2006-08.  (They  are  available  in  pdf  format  at  http://eccentric.mae.comell.edu/~tcg/pubs.html.) 

[1]  B.  Merci,  D.  Roekaerts,  B.  Naud  and  S.B.  Pope  (2006)  “Comparative  study  of  micro¬ 
mixing  models  in  transported  scalar  PDF  simulations  of  turbulent  non-premixed  bluff  body 
flames,”  Combustion  and  Flame,  146,  109-130. 

[2]  R.L.  Gordon,  A.R.  Masri,  S.B.  Pope  and  G.M.  Goldin  (2007a)  “A  Numerical  Study  of 
Auto-ignition  in  Turbulent  Lifted  Flames  Issuing  into  a  Vitiated  Co-flow,”  Combustion 
Theory  and  Modelling.  1 1  (3),  351-376. 

[3]  R.L.  Gordon,  A.R.  Masri,  S.B.  Pope  and  G.M.  Goldin  (2007b)  “Transport  budgets  in 
turbulent  lifted  flames  of  methane  auto-igniting  in  a  vitiated  co-flow,”  Combustion  and 
Flame,  151  (3),  495-511. 

[4]  M.R.H.  Sheikhi,  P.  Givi  and  S.B.  Pope  (2007)  “Velocity-scalar  filtered  mass  density 
function  for  large  eddy  simulation  of  turbulent  reacting  flows,”  Physics  of  Fluids  1 9  (9), 
095106. 

[5]  H.  Wang  and  S.B.  Pope  (2008a)  “Lagrangian  investigation  of  local  extinction,  re-ignition 
and  auto-ignition  in  turbulent  flames,”  Combustion  Theory  and  Modelling,  Combustion 
Theory  and  Modelling  12(5),  857-882. 

[6]  H.  Wang  and  S.B.  Pope  (2008b)  “Time  averaging  strategies  in  the  finite-volume/particle 
algorithm  for  the  joint  PDF  equation  of  turbulent  reactive  flows,”  Combustion  Theory  and 
Modelling,  Combustion  Theory  and  Modelling  12(3),  529-544. 

[7]  L.  Lu,  S.  R.  Lantz,  Z.  Ren  and  S.  B.  Pope  (2008)  “Computationally  efficient 
implementation  of  combustion  chemistry  in  parallel  PDF  calculations,”  Journal  of 
Computational  Physics,  (submitted). 

[8]  P.P.  Popov,  R.  McDermott  and  S.B.  Pope  (2008)  “An  accurate  Time  Advancement 
Algorithm  for  Particle  Tracking,”  Journal  of  Computational  Physics,  Volume  227,  Issue  20, 
pp.  8792-8806. 

[9]  H.  Wang,  P.  P.  Popov,  and  S.  B.  Pope  (2008)  “Weak  Second  Order  Splitting  Schemes  for 
Lagrangian  Monte  Carlo  Particle  Methods  for  the  Composition  PDF/FDF  Transport 
Equations,”  Journal  of  Computational  Physics  ,  (submitted). 

OTHER  REFERENCES 

[1]  R.S.  Barlow  and  J.H.  Frank  (1998)  “Effects  of  Turbulence  on  Species  Mass  Fractions  in 
Methane/Air  Jet  Flames,”  Proc.  Combust.  Inst.  27,  1087-1095. 

[2]  J.  Xu,  and  S.B.  Pope  (2000)  “PDF  calculations  of  turbulent  nonpremixed  flames  with  local 
extinction,”  Combust.  Flame  123,  281-307. 

[3]  Q.  Tang,  J.  Xu  and  S.B.  Pope  (2000)  “PDF  calculations  of  local  extinction  and  NO 
production  in  piloted-jet  turbulent  methane/air  flames,”  Proceedings  of  the  Combustion 
Institute,  28,  133-139. 


16 


[4]  R.P.  Lindstedt,  S.A.  Louloudi,  and  E.M.  Vaos  (2000)  “Joint  scalar  probability  density 
function  modeling  of  pollutant  formation  in  piloted  turbulent  jet  diffusion  flames  with 
comprehensive  chemistry,”  Proc.  Combust.  Inst.  28,  149-156. 

[5]  R.S.  Barlow,  A.N.  Karpetis,  J.H.  Frank,  J.Y.  Chen,  (2001)  “Scalar  profiles  and  NO 
formation  in  laminar  opposed-flow  partially  premixed  methane/air  flames,”  Combust. 
Flame  127,  2102-2118. 

[6]  R.  Cabra,  T.  Myhrvold,  J.Y.  Chen,  R.W.  Dibble,  A.N.  Karpetis,  and  R.S.  Barlow  (2002) 
“Simultaneous  laser  raman-rayleigh-lif  measurements  and  numerical  modeling  results  of  a 
lifted  turbulent  H2/N2  jet  flame  in  a  vitiated  coflow,”  Proc.  Combust.  Inst.  29,  1881-1888. 

[7]  A.R.  Masri,  R.  Cao,  S.B.  Pope  and  G.M.  Goldin  (2004)  “PDF  Calculations  of  Turbulent 
Lifted  Flames  of  H2/N2  issuing  into  a  vitiated  co-flow,”  Combust.  Theory  Modelling,  8,  1- 
22. 

[8]  R.  Cabra,  J.-Y.  Chen,  R.W.  Dibble,  A.N.  Karpetis  and  R.S.  Barlow  (2005)  “Lifted 
methane-air  jet  flames  in  a  vitiated  coflow,”  Combustion  and  Flame,  143,  491-506. 

[9]  R.  Cao  and  S.B.  Pope  (2005)  “The  influence  of  chemical  mechanisms  on  PDF  calculations 
of  nonpremixed  piloted  jet  flames,”  Combustion  and  Flame,  143,  450-470. 

[10]  R.  Cao,  H.  Wang  and  S.B.  Pope  (2005)  “The  effect  of  mixing  models  in  PDF  calculations 
of  piloted  jet  flames,”  Proc.  Combust.  Inst.,  31,  1543-1550. 

[10]  R.  Cao,  S.B.  Pope  and  A.R.  Masri  (2005)  “Turbulent  lifted  flames  in  a  vitiated  coflow 
investigated  using  joint  PDF  calculations,”  Combustion  and  Flame,  142,  549-568. 

[11]  S.B.  Pope  (1997)  “Computationally  efficient  implementation  of  combustion  chemistry 
using  in  situ  adaptive  tabulation,”  Combust.  Theory  Modelling  1, 41-63. 

[11]  L.  Lu,  Z.  Ren,  S.R.  Lantz,  V.  Raman,  S.B.  Pope  and  H.  Pitsch  (2005)  “Investigation  of 
strategies  for  the  parallel  implementation  of  ISAT  in  LES/FDF/1SAT  computations,”  4th 
Joint  Meeting  of  the  U.  S.  Sections  of  the  Combustion  Institute,  Philadelphia,  PA. 

[12]  Z.  Ren,  and  S.B.  Pope  (2004)  “An  investigation  of  the  performance  of  turbulent  mixing 
models,”  Combust.  Flame  136,  208-216. 

[13]  FI.  Wang,  S.B.  Pope  (2007)  “Comparison  of  detailed  reaction  mechanism  in  non-premixed 
combustion,”  5th  Joint  Meeting  of  the  U.  S.  Sections  of  the  Combustion  Institute,  UCSD, 
San  Diego,  March  26-28. 

[14]  R.  McDermott  and  S.B.  Pope  (2008),  “The  Parabolic  Edge  Reconstruction  Method 
(PERM)  for  Lagrangian  Particle  Advection,”  Journal  of  Computational  Physics,  227  (10), 
5447-5491. 

[15]  R.  Cao  and  S.B.  Pope  (2003)  “Numerical  Integration  of  Stochastic  Differential  Equations: 
Weak  Second-Order  Mid-Point  Scheme  for  Applications  in  the  Composition  PDF 
Method,”  Journal  of  Computational  Physics,  185,  194-212. 

[16]  A.  Tocino  and  J.  Vigo-Aguiar  (2002)  “Weak  Second  Order  Conditions  for  Stochastic 
Runge-Kutta  Methods,”  SIAM  Journal  on  Scientific  Computing,  24  (2),  507-523. 

[17]  P.P.  Popov  and  S.B.  Pope  (2008)  “Stochastic  Particle  Advection  in  Hybrid  Large  Eddy 
Simulation/Filtered  Density  Function  Methods,”  61st  Annual  Meeting  of  the  American 
Physical  Society,  Division  of  Fluid  Dynamics,  San  Antonio,  November  23-25,  2008. 

[18]  S.B.  Pope  (1985)  “PDF  methods  for  turbulent  reactive  flows,”  Progress  in  Energy  and 
Combustion  Science,  11,  119-192. 

[19]  P.E.  Kloeden,  E.  Platen  (1992)  “Numerical  Solution  of  Stochastic  Differential  Equations,” 
Springer-Verlag,  Berlin. 

[20]  P.J.  Roache  (2002)  “Code  verification  by  the  method  of  manufactured  solutions,”  J.  Fluids 
Eng.  124(1),  4-10. 

[21]  S.  Subramaniam  and  S.B.  Pope  (1998)  “A  mixing  model  for  turbulent  reactive  flows  based 
on  Euclidean  minimum  spanning  trees,”  Combustion  and  Flame,  1 15,  487--514. 


17 


PERSONNEL  SUPPORTED 


Prof.  Stephen  B.  Pope,  PI 
Liuyan  Lu,  Ph.  D.  student 
Haifeng  Wang,  Ph.D.  student 
Pavel  Popov,  Ph.D.  student 

DEGREES  GRANTED 

Liuyan  Lu,  Ph.D.,  Mechanical  Engineering,  2008. 

MEETING  PARTICIPATION 

11th  International  Conference  on  Numerical  Combustion,  Granada,  Spain,  April  23-26,  2006. 
Two  contributed  talks. 

27th  Annual  Combustion  Research  Meeting,  DOE/BES,  Wintergreen,  VA,  May  30-June  2,  2006. 

ARO/AFOSR  Contractors’  Meeting,  in  Chemical  Propulsion,  Arlington  ,  VA,  June  12-14,  2006. 

7th  World  Congress  on  Computational  Mechanics,  Los  Angeles,  CA,  July  16-22,  2006.  Invited 
keynote  lecture. 

8th  International  Workshop  on  Measurement  and  Computation  of  Turbulent  Nonpremixed 
Flames,  Heidelberg,  Germany,  August  3-5,  2006. 

31v  International  Symposium  on  Combustion,  Heidelberg,  Germany,  August  6-11,  two 
contributed  papers. 

59th  Annual  Meeting  of  the  American  Physical  Society,  Division  of  Fluid  Dynamics,  November 
19-21,  2006,  Tampa,  FL.  Five  contributed  talks. 

Johns  Hopkins  University,  April  13,  2007,  Seminar. 

NSF  Panel,  Washington,  D.C.,  May  1, 2007. 

28th  Annual  Combustion  Research  Meeting,  DOE/BES,  Airlie,  VA,  May  29-June  1, 2007. 
ARO/AFOSR  Contractors’  Meeting  in  Chemical  Propulsion,  Boulder,  CO,  June  11-13,  2007. 
Royal  Society,  New  Fellows  Seminar,  London,  July  11-12,  2007. 

ECCOMAS  Conference  on  Computational  Combustion,  Delft,  Netherlands,  July  18-20,  2007. 
Invited  talk. 

Air  Force  Research  Laboratories,  Wright-Patterson  Air  Force  Base,  Dayton,  OH,  August  8, 
2007.  Seminar  and  visit. 


18 


Pennsylvania  State  University,  Air  Products  Distinguished  lecture,  October  9,  2007. 

APS/DFD  Annual  Meeting,  Salt  Lake  City,  November  18-20,  2007. 

Delft  University,  Netherlands,  The  Burgers  Lecture,  January  10,  2008. 

SIAM  12'  International  Conference  on  Numerical  Combustion,  Monterey,  March  31 -April  2, 
2008. 

DOE  29th  Annual  Combustion  Research  Conference,  Airlie,  May  27-30,  2008. 

Symposium  on  Fluid  Science  and  Turbulence,  Johns  Hopkins  University,  May  30-31, 2008. 

AFOSR-ARO  Contractors  Meeting:  Chemical  Propulsion,  Arlington,  June  16-18,  2008. 

TNF9:  Ninth  International  Workshop  on  Measurement  and  Computation  of  Turbulent 
Nonpremixed  Flames,  Montreal,  July  31 -August  2,  2008. 

32nri  International  Symposium  on  Combustion,  Montreal,  August  3-8,  2008. 

DLES7:  ERCOFTAC  Workshop  on  Direct  and  Large-Eddy  Simulation,  Trieste,  Italy,  September 
8-10,  2008. 

DNS  and  LES  of  Reacting  Flows,  Maastricht,  Netherlands,  October  22-24,  2008. 

61st  Annual  Meeting  of  the  American  Physical  Society,  Division  of  Fluid  Dynamics,  San 
Antonio,  November  23-25,  2008. 

INTERACTION  WITH  AFRL 

At  the  ARO/AFOSR  Contractors’  Meeting  (June  12-14,  2006)  I  discussed  with  Dr.  J.  T.  Edwards 
the  possible  application  to  AFRL  scramjet  computations  of  chemistry  dimension  reduction 
techniques  developed  at  Cornell.  Following  Dr.  Edwards’  suggestion,  I  then  had  an  e-mail 
correspondence  with  Dr.  Campbell  Carter  to  provide  further  information  on  the  technique,  and  to 
suggest  how  it  might  be  applied.  AFRL  are  currently  considering  the  possibilities. 

On  August  8,  2007,  the  PI  visited  AFRL  and  presented  a  seminar  on  “Computational  Modeling 
of  Turbulence-Chemistry  Interactions  in  Combustion.”  In  addition  to  touring  several  of  the 
facilities,  he  discussed  ongoing  research  with  Cam  Carter,  Mel  Roquemore,  Doug  Davis,  Susan 
Cox-Stouffer,  among  others.  This  included  discussions  on  the  modeling  of  turbulence-chemistry 
interactions  in  computations  of  scramjet  engines.  Since  the  AFOSR-ARO  Contractors  meeting 
in  June  2008,  the  PI  has  contacted  Drs.  Campbell  Carter,  Bish  Ganguly,  Skip  Williams  and 
James  Gord  in  order  to  continue  the  dialog  on  combustion  modeling  at  AFRL. 

On  October  16,  2008,  the  PI  hosted  a  visit  to  Cornell  by  Dr.  Barry  Kiel  of  AFRL.  Dr.  Kiel  made 
a  presentation  on  augmentor  research,  and  had  discussions  with  the  PI  and  his  research  group  on 
this  topic.  Subsequently,  the  PI  provided  Dr.  Kiel  with  a  written  assessment  of  the  role  of 
LES/FDF  methods  in  augmentor  simulations. 


19 


TECHNOLOGY  TRANSITIONS  AND  TRANSFERS 


Over  the  years,  the  PI  has  successfully  transferred  several  aspects  of  PDF  methods  and  their 
numerical  implementation  to  Fluent,  Inc.  As  Fluent  is  the  most  widely  used  CDF  code,  this  is  an 
effective  way  to  transition  PDF  capabilities  to  the  aerospace  (and  other)  industries.  For  example, 
the  combustor  design  group  at  GE  Aircraft  Engines  is  using  the  PDF  capability  in  Fluent.  There 
is,  however,  no  direct  technology  transfer  to  document  this  period. 

HONORS  AND  AWARDS 

Zeldovich  Gold  Medal,  The  Combustion  Institute 
Fellow  of  the  Royal  Society 

Fellow  of  the  American  Academy  of  Arts  and  Sciences 
Fellow  of  the  American  Physical  Society 
Fellow  of  the  Institute  of  Physics 
Associate  Fellow  of  A1AA 

INVENTIONS  AND  PATENTS 

There  were  no  inventions  or  patents  during  the  reporting  period. 


20 


Combustion  Theory  and  Modelling 
Vol.  12,  No.  5,  October  2008,  857-882 


©Taylor  &  Francis 

Tartof4Hand»Cf««i 


Lagrangian  investigation  of  local  extinction,  re-ignition  and 
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Lagrangian  PDF  investigations  are  performed  of  the  Sandia  piloted  flame  E  and  the  Cabra 
H2/N2  lifted  flame  to  help  develop  a  deeper  understanding  of  local  extinction,  re-ignition  and 
auto-ignition  in  these  flames,  and  of  the  PDF  models’  abilities  to  represent  these  phenomena. 
Lagrangian  particle  time  series  are  extracted  from  the  PDF  model  calculations  and  are  analyzed. 

In  the  analysis  of  the  results  for  flame  E,  the  particle  trajectories  are  divided  into  two  groups: 
continuous  burning  and  local  extinction.  For  each  group,  the  trajectories  are  further  sub-divided 
based  on  the  particles’  origin:  the  fuel  stream,  the  oxidizer  stream,  the  pilot  stream,  and  the 
intermediate  region.  The  PDF  calculations  are  performed  using  each  of  three  commonly  used 
models  of  molecular  mixing,  namely  the  EMST,  1EM  and  modified  Curl  mixing  models. 

The  calculations  with  different  mixing  models  reproduce  the  local  extinction  and  re-ignition 
processes  observed  in  flame  E  reasonably  well.  The  particle  behavior  produced  by  the  1EM  and 
modified  Curl  models  is  different  from  that  produced  by  the  EMST  model,  i.e.,  the  temperature 
drops  prior  to  (and  sometimes  during)  re-ignition.  Two  different  re-ignition  mechanisms  arc 
identified  for  flame  E:  auto-ignition  and  mixing-reaction.  In  the  Cabra  H2/N2  lifted  flame,  the 
particle  trajectories  are  divided  into  different  categories  based  on  the  particles’  origin:  the  fuel 
stream,  the  oxidizer  stream,  and  the  intermediate  region.  The  calculations  reproduce  the  whole 
auto-ignition  process  reasonably  well  for  the  Cabra  flame  Four  stages  of  combustion  in  the 
Cabra  flame  are  identified  in  the  calculations  by  the  different  mixing  models,  i  e.,  pure  mixing, 
auto-ignition,  mixing-ignition,  and  fully  burnt,  although  the  individual  particle  behavior  by 
the  IEM  and  modified  Curl  models  is  different  from  that  by  the  EMST  model.  The  relative 
importance  of  mixing  and  reaction  during  re-ignition  and  auto-ignition  are  quantified  for  the 
IEM  model. 

Key  words:  PDF  method,  particle  tracking,  particle  trajectories,  mixing  models,  local  extinc¬ 
tion,  re-ignition,  auto-ignition 


1.  Introduction 

Flame  extinction  and  ignition  are  fundamental  phenomena  in  combustion  problems.  The  occur¬ 
rence  of  extinction  and  ignition  in  turbulent  reactive  flows,  due  to  intensive  non-linear  turbulence- 
chemistry  interactions,  is  a  challenge  to  modem  turbulent  combustion  models.  The  probability 
density  function  (PDF)  transport  equation  method  [1-3]  is  increasingly  found  to  be  able  to  ac¬ 
count  accurately  for  the  turbulence- chemistry  interactions,  e.g.,  local  extinction  and  re-ignition 
[4,  5].  In  engineering  practice,  the  statistics  of  the  turbulent  velocity  and  composition  fields  are 
of  primary  concern  in  the  context  of  Reynolds  averaged  Navier-Stokes  (RANS)  simulations.  In 
PDF  methods,  the  modeled  PDF  transport  equation  is  usually  solved  numerically  by  a  Lagrangian 
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particle  method,  and  much  more  information  can  be  extracted  from  the  particle  properties.  Par¬ 
ticle  scatter  plots  (or  joint  PDFs)  contain  the  most  detailed  information  about  the  distribution  of 
properties  at  a  given  position  and  time.  All  previous  PDF  calculations  of  turbulent  flames  have 
focused  on  Eulerian  statistics  and  their  comparison  with  experimental  data,  e.g.,  conditional  or 
unconditional  statistics  of  compositions,  particle  scatter  plots  and  conditional  PDFs  of  compo¬ 
sitions.  Scatter  plots  of  particle  properties  are  able  to  illustrate  qualitatively  different  kinds  of 
complicated  turbulent  combustion  phenomena  such  as  local  extinction  and  re-ignition.  To  ex¬ 
plore  the  PDF  calculation  results  comprehensively,  and  to  understand  the  turbulence-chemistry 
interactions  more  deeply,  here  we  extract  and  analyze  Lagrangian  time  series  of  particle  proper¬ 
ties  from  the  PDF  calculations.  The  particle  trajectories  are  presented  to  illustrate  the  dynamic 
evolution  of  complicated  turbulent  combustion  processes,  i.e.,  local  extinction  and  re-ignition 
in  the  turbulent  non-premixed  piloted  jet  flame,  and  auto-ignition  in  the  turbulent  lifted  jet 
flame. 

Lagrangian  properties  are  important  physical  properties  relevant  not  only  to  the  PDF  method, 
but  also  to  real  turbulence  and  combustion  problems.  Many  Lagrangian  investigations  have 
been  performed  of  turbulence  using  direct  numerical  simulations  (DNS).  Yeung  [6,  7]  studied 
the  Lagrangian  characteristics  of  turbulence  and  passive  scalar  transport  in  stationary  isotropic 
turbulence  with  uniform  mean  scalar  gradients.  The  Lagrangian  properties  of  the  scalars  investi¬ 
gated  are  important  to  molecular  mixing  models.  Mitarai  et  al.  [8]  performed  DNS  of  an  idealized 
non-premixed  flame  in  decaying  isotropic  turbulence  for  conditions  w'here  flame  extinction  and 
re-ignition  occur.  In  that  work,  the  fluid  particles  are  tracked  to  investigate  flame  extinction  and 
re-ignition.  Different  categories  of  particles  are  identified,  e.g.,  continuous  burning  and  local 
extinction.  Also  investigated  are  Lagrangian  properties  of  the  conditional  scalar  diffusion,  which 
appears  as  an  unclosed  term  in  the  PDF  transport  equation.  The  same  methodology  is  used  to  in¬ 
vestigate  the  performance  of  flamelet  models  [9]  and  the  performance  of  different  mixing  models 
[10].  Sripakagonmtf/.  [1 1]  performed  Lagrangian  flame  element  tracking  along  the  stoichiomet¬ 
ric  surface  in  decaying  isotropic  turbulence  to  investigate  flame  extinction  and  re-ignition.  Three 
major  scenarios  of  re-ignition  in  non-premixed  combustion  are  identified,  i.e.,  the  independent 
flamelet  scenario,  re-ignition  via  edge  flame  propagation,  and  re-ignition  through  engulfment 
by  hot  neighbouring  fluid.  These  Lagrangian  investigations  are  helpful  to  provide  insights  into 
the  dynamic  evolution  of  turbulent  combustion  processes.  Because  of  the  formidable  practical 
difficulties  there  are  no  experimental  data  on  Lagrangian  quantities  in  turbulent  reactive  flows. 
Experimental  data  are,  however,  becoming  available  on  Lagrangian  velocity  and  acceleration 
statistics  in  non- reactive  flows  [12-14]. 

The  Lagrangian  PDF  method  [2]  represents  the  turbulent  flow,  transport  and  reaction  processes 
via  the  time  evolution  of  nominal  Monte  Carlo  particles  representing  the  joint  PDF  of  velocity, 
turbulence  frequency  and  compositions.  Eulerian  statistics  obtained  from  PDF  calculations  have 
been  explored  extensively  before,  and  are  generally  found  to  be  in  good  agreement  with  the 
experimental  data.  It  would  be  valuable  to  extract  the  Lagrangian  time  series  which  contains 
the  whole  history  of  the  particle  evolution  in  the  multi-dimensional  sample  space.  This  work  is 
dedicated,  as  a  first  effort,  to  explore  the  Lagrangian  properties  of  the  Monte  Carlo  particles  in  the 
PDF  simulations  of  turbulent  flames  containing  local  extinction  and  re-ignition,  and  auto-ignition 
The  two  flames  studied  are  the  Sandia  non-premixed  piloted  jet  flame  E  [  1 5],  and  the  lifted  H2/N2 
jet  flame  in  a  vitiated  coflow  [16],  referred  to  as  the  Cabra  flame. 

In  PDF  methods,  the  closed  form  of  the  chemical  reaction  source  term  facilitates  the  exact 
treatment  of  detailed  combustion  chemistry.  The  modeling  of  the  unclosed  molecular  mixing 
term  in  the  PDF  equation  remains  one  of  the  major  efforts  of  model  development.  Three  mixing 
models  are  extensively  used,  i.e.  the  Euclidean  minimum  spanning  tree  (EMST)  model  [21],  the 
interaction  by  exchange  with  the  mean  (IEM)  model  [18]  (or  the  least-mean-square  estimator 
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(LMSE)  model  [19]),  and  the  modified  Curl  model  [20,  22].  All  the  mixing  models  can  represent 
local  extinction  and  re-ignition  to  some  extent  [4,  5,  23],  although  the  EMST  model  is  usually 
thought  to  be  superior  in  this  respect.  The  simplicity  of  the  1EM  and  modified  Curl  models  makes 
them  quite  popular.  In  spite  of  the  complexity  of  the  EMST  model,  the  public  availability  of  a 
FORTRAN  implementation  [24]  makes  it  easy  to  use.  A  desirable  property  of  mixing  models 
is  “localness”  [21].  All  the  models  are  local  in  physical  space;  only  the  EMST  model  is  local 
in  composition  space;  and  none  is  local  in  velocity  space.  There  is  some  recent  progress  in 
the  development  of  more  sophisticated  mixing  models,  e.g.,  the  multiple  mapping  conditioning 
(MMC)  model  [25], and  the  interaction  by  exchange  with  the  conditional  mean  (LECM)  mixing 
model  [26-29],  which  is  local  in  velocity  space.  The  present  work  focuses  on  the  three  traditional 
mixing  models  (  EMST,  IEM  and  modified  Curl)  and  evaluates  their  relative  performance  from 
the  Lagrangian  viewpoint. 

The  primary  aim  of  PDF  methods  is  to  calculate  accurately  one-point,  one-time  Eulerian 
quantities.  It  is  well  understood  [  1  ]  that  this  objective  can  be  achieved  using  stochastic  Lagrangian 
models,  even  if  the  multi-time  properties  of  the  models  are  not  physically  accurate.  While  the 
multi-time  behavior  of  the  stochastic  models  for  position  and  velocity  are  physically  realistic, 
those  of  the  mixing  models  are  not.  For  example  Curl’s  model  involves  jumps  in  compositions; 
and  (for  the  simplest  homogeneous  turbulence)  the  IEM  model  yields  a  deterministic  relaxation  to 
the  mean,  with  no  fluctuations  along  Lagrangian  trajectories.  Hence,  while  this  study  is  valuable 
in  shedding  light  on  the  models’  behavior  and  performance,  a  close  correspondence  between  the 
models’  Lagrangian  trajectories  and  those  in  the  flame  (could  they  be  measured!)  should  not  be 
expected. 

The  remaining  sections  of  this  paper  are  organized  as  follows.  In  Section  2,  Eulerian  scatter 
plots  of  particle  properties  are  presented  for  the  Sandia  piloted  flame  E  [15]  and  for  the  Cabra 
H2/N2  lifted  jet  flame  [16].  The  illustrations  of  local  extinction  and  re-ignition,  and  of  auto- ignition 
are  reviewed  by  reference  to  the  Eulerian  particle  data.  The  limitations  of  the  Eulerian  data  are 
discussed.  In  Section  3,  the  particle  tracking  and  particle  sampling  procedures  are  presented.  In 
Sections  4  and  5  (for  the  Sandia  flame  E  and  the  Cabra  lifted  flame,  respectively),  the  Lagrangian 
time  series  obtained  from  the  PDF  calculation  using  the  different  mixing  models  are  analyzed  to 
study  the  models’  representation  of  extinction,  re-ignition  and  auto-ignition.  The  relative  roles 
of  mixing  and  reaction  during  re-ignition  and  auto-ignition  are  quantified  for  the  IEM  model  in 
Section  6,  Conclusions  are  drawn  in  the  final  section. 

2.  Particle  calculations  and  Eulerian  scatter  plots 

Comprehensive  PDF  model  investigations  of  the  Sandia  piloted  flames  and  the  Cabra  H2/N2  lifted 
jet  flame  have  been  described  elsewhere  [4,  5,  16,  23,  30-32].  In  this  section,  PDF  calculations 
of  the  Sandia  flame  E  and  the  Cabra  lifted  jet  flame  are  repeated  to  review  the  Eulerian  particle 
scatter  plots.  As  in  [23,  30],  a  PDF  code  called  HYB2D  is  used,  in  which  a  hybrid  finite  volume 
(FV)/particle  algorithm  is  implemented  for  solving  the  joint  PDF  transport  equation  of  the 
velocity,  turbulence  frequency  and  compositions  [34].  The  details  of  the  simulations  for  the  Sandia 
flame  E  and  the  Cabra  lifted  jet  flame  are  identical  to  those  in  [35]  and  in  [31],  respectively,  and 
are  simply  summarized  in  Table  1 .  (Quantities  not  listed  in  Table  1  can  be  found  in  [3 1  ]  and  [35].) 
Different  values  of  the  mixing  model  constant  C<p  are  specified  for  the  different  mixing  models  in 
the  calculation  of  the  Sandia  flame  E  in  order  to  achieve  a  stable  burning  flame  with  roughly  the 
same  amount  of  local  extinction  [23,  35]  as  observed  experimentally.  Similarly,  different  coflow 
temperatures  are  used  in  the  calculation  of  the  Cabra  lifted  flame  with  different  mixing  models 
in  order  to  produce  approximately  the  same  flame  lift-off  height  as  observed  experimentally 
[31]. 


860 


H.  Wang  and  S.B.  Pope 


Table  1  Details  of  the  simulations  for  the  Sandia  flame  E  and  the  Cabra  lifted  H2/N2  jet  flame 


Model  parameters 

Sandia  flame  E 

Cabra  lifted  flame 

Turbulence  frequency  model 
constant  [17],  C^i 

0.7 

0.65 

Number  of  particles  per  cell,  Npc 

100 

100 

Chemistry 

GRl-Mech  3.0  [37] 

H2-O2  mechanism  [38] 

1SAT  error  tolerance  [36],  etoi 

5.0  x  10"5 

6.25  x  10‘6 

Grid  size 

96  x  96 

96  x  96 

EMST  1EM  modified  Curl 

EMST  1EM 

modified  Curl 

Mixing  model  constant,  Q> 

1.5  2.7  3.3 

1.5  1.5 

1.5 

Coflow  temperature  in  Cabra 
lifted  flame,  Tc( K) 

1033  1036 

1036 

Conventionally,  the  output  from  PDF  calculations  is  Eulenan  data  for  analysis  at  a  fixed 
time  (when  the  statistically  stationary  state  has  been  reached)  and  at  different  locations.  The 
conditional  and  unconditional  statistics  of  the  Eulerian  data  have  been  discussed  extensively 
elsewhere  [4,  5,  23,  30,  31]  and  will  not  be  repeated  here.  In  this  section,  we  review  the  scat¬ 
ter  plots  of  the  particle  temperature  versus  the  mixture  fraction  in  the  Sandia  flame  E  and  the 
Cabra  lifted  flame.  As  previously  discussed  [4],  it  is  difficult  to  make  a  rigorous  quantitative 
comparison  between  scatter  plots  from  experiment  and  model  calculations,  because  of  differ¬ 
ences  in  the  sampling  and  weighting  of  particles.  Nevertheless,  this  comparison  is  useful  in 
assessing,  at  least  qualitatively,  the  ability  of  the  models  to  represent  the  phenomena  observed 
experimentally. 

Figure  1  shows  the  scatter  plots  of  the  particle  temperature  against  the  mixture  fraction  at 
different  axial  locations  in  the  Sandia  flame  E  from  measurements  [15],  and  from  the  PDF 
simulations  with  different  mixing  models.  The  axial  distance  x  is  shown  as  x /D,  where  D  is  the 
diameter  of  the  fuel  jet.  Two  laminar  flame  temperature  profiles  (dashed  lines  in  Figure  1)  are 
also  shown  for  reference.  The  laminar  calculations  are  conducted  by  using  OPPDIF  [39]  with  two 
strain  rates,  a  =  10  s~]  and  310s-1,  and  in  the  latter  case,  the  temperature  profile  is  shifted  down 
by  300  K.  Following  [8],  we  use  this  shifted  temperature  profiles  with  a  =310  s_1  as  a  simple 
criterion  to  distinguish  between  burning  particles  (above  the  line)  and  extinguished  particles 
(below  the  line).  For  simplicity,  we  call  this  line  the  ‘‘extinction  line”,  and  the  region  above  the 
“burning  region”,  and  the  region  below  the  “extinction  region”.  There  is  of  course  an  extinction 
limit  of  strain  rate  ae  ( ae  ^  376  s-1  for  the  current  case)  and  its  corresponding  temperature  profile 
in  steady  opposed  laminar  non-premixed  jet  flames,  but  we  prefer  not  to  use  this  extinction  limit 
as  our  criterion.  The  extinction  limit  af  applies  only  to  steady  laminar  flames.  In  the  unsteady 
case  (e.g.,  laminar  flames  subject  to  the  oscillation  of  strain  rate),  however,  the  instantaneous 
strain  rate  can  exceed  the  extinction  limit,  without  the  flame  being  extinguished  [40,  41].  Using 
the  extinction  limit  of  the  steady  laminar  flame  will  somewhat  over-estimate  the  amount  of  local 
extinction  in  this  turbulent  flame.  Flence,  as  in  [8],  we  use  the  shifted  temperature  profiles  as  a 
more  conservative  extinction  limit.  This  criterion  is  somewhat  arbitrary,  but  it  is  helpful  for  the 
qualitative  analysis  of  the  local  extinction  and  re-ignition  reported  below. 

From  the  experimental  scatter  plots  in  Figure  1,  it  may  be  seen  that  the  number  of  the 
particles  below  the  extinction  lines  decreases  with  increasing  the  axial  distance  from  x/D  =  7.5 
to  45,  indicating  the  evolution  from  local  extinction  to  re-ignition.  The  PDF  calculations  with 
the  three  mixing  models  qualitatively  reproduce  this  process  to  some  extent.  The  scatter  plots 
of  the  Eulerian  particle  data  visually  illustrate  the  level  of  local  extinction  at  different  locations. 
However,  the  Lagrangian  evolution  of  the  particle  properties  is  not  evident.  Where  do  the  locally 
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Figure  1.  Scatter  plots  of  temperature  against  the  mixture  fraction  at  the  axial  locations  x/D  =  7.5,  15, 
30,  and  45  in  the  Sandia  piloted  flame  E  from  experimental  data  and  from  PDF  calculations  using  three 
different  mixing  models  (Open  circles:  conditional  mean  temperature;  Dashed  lines:  temperature  profiles 
in  the  opposed-jet  laminar  flame  with  strain  rate  a  =  10s-1  (upper  lines);  and  with  strain  rate  a  —  3 10s_l 
(lower  lines),  shifted  down  by  300  K.) 


extinguished  particles  come  from?  How  do  they  return  to  the  burning  region  (in  composition 
space)?  How  do  different  mixing  models  cause  the  particles  to  move  in  the  composition  space? 
The  current  Eulerian  data  cannot  answer  these  questions. 

We  now  turn  our  attention  to  the  Cabra  H2/N2  lifted  jet  flame.  Figure  2  shows  the  scatter  plots 
of  temperature  versus  mixture  fraction  at  different  axial  locations  in  the  flame.  The  equilibrium 
state  calculated  by  using  EQUIL  [42]  is  also  shown  in  the  plots  for  reference.  The  initial  enthalpy 
h  and  the  species  mass  fractions  Y  for  the  equilibrium  calculation  are  taken  to  be  linear  in  the 
mixture  fraction  £  space,  i.e. 


h(Z)  =  h0\  -  (h0x  -  Afii)  •  £ 


(1) 


Y(£)  =  Yox  —  (Yox  —  Yfu)  ♦  £ 


(2) 


where  the  subscript  “fu”  and  “ox”  denote  the  fuel  stream  and  the  oxidizer  stream,  respectively. 
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Figure  2.  Scatter  plots  of  particle  temperature  against  mixture  fraction  at  the  axial  locations  x /  D  =9,  1 1, 
14,  and  26  in  the  Cabra  H2/N2  lifted  flame  from  experimental  data  and  from  PDF  calculations  using  three 
different  mixing  models  (Open  circles:  conditional  mean  temperature;  Dashed  lines:  equilibrium  state.) 


From  the  experimental  data  at  x/D  =  9  shown  in  Figure  2,  it  may  be  seen  that  the  particles 
lie  dominantly  on  the  mixing  line,  with  just  a  few  rare  particles  with  higher  temperature.  This 
combustion  stage  (x/D  ^  9)  is  called  pure  mixing.  For  x/D  ^  1 1,  the  particles  leave  the  mixing 
line  gradually,  indicating  an  ignition  process.  By  x/D  =  26,  almost  all  the  particles  have  reached 
the  fully  burnt  state,  close  to  the  equilibrium  line.  The  PDF  calculations  using  the  three  mixing 
models  predict  the  mixing-ignition  processes  reasonably  well.  However,  the  scatter  plots  of 
temperature  are  different  for  the  different  mixing  models.  Similar  questions  arise,  e.g.,  how  do 
the  different  mixing  models  cause  the  particles  to  evolve  through  the  mixing  stage  to  the  burning 
state?  The  Eulerian  data  cannot  answer  such  questions. 


3.  Lagrangian  particle  tracking 

The  limitations  of  the  Eulerian  data  from  PDF  calculations  are  evident.  In  this  section,  we  discuss 
the  extraction  of  Lagrangian  time  series  from  the  PDF  calculations,  and  these  are  analyzed  for 
the  Sandia  flame  E  and  the  Cabra  lifted  flame  in  the  following  two  sections. 

In  the  Lagrangian  PDF  method  [2],  the  modeled  transport  equation  for  the  joint  PDF  of  the 
velocity  U,  turbulence  frequency  co ,  and  the  compositions  <j>  is  solved  by  a  Monte  Carlo  particle 
method.  A  large  number  of  Monte  Carlo  particles  are  released  into  the  computational  domain 
according  to  the  Eulerian  PDF  initially.  Each  particle  carries  a  full  set  of  the  fluid  properties,  i.e.. 
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U*,  a mass  m*,  locations  x*  and  4> *  etc.  The  evolution  of  the  joint  PDF  is  represented  by  the 
movement  of  the  particles  in  the  multi-dimensional  space  governed  by  the  following  stochastic 


differential  equations  [1-3] 

dx*  =  U*dt ,  (3) 

du;  =  -  Q  +  ^Co)  -  u,)dt  +  (C0kn)'/2dW,,  (4) 

dw*  =  -<S)dt  -  SuQco'dt  +  {2Cw3C<04wQ(o*)'/2  dW,  (5) 

^  =Ma(t)  +  Sa(<l>*(t)),  (6) 

dt 


where  p  and  p  are  the  fluid  density  and  pressure,  respectively;  “( )”  denotes  the  conventional  mean; 

denotes  the  Favre  mean;  Co,  3  and  C ^  are  model  constants;  W  is  an  isotropic  vector- valued 
Wiener  process;  W  is  another  independent  Wiener  process;  k  is  the  turbulent  kinetic  energy;  S(l) 
and  Sa  are  the  source  term  for  cd  and  the  reaction  source  term  for  <pa ,  respectively;  MaU)  denotes 
the  mixing  model;  Q  is  the  conditional  mean  turbulence  frequency  defined  as 


Q  =  Cq 


(p*co*\co*  ^  co) 

Tp) 


(7) 


where  the  constant  Cq  is  chosen  so  that  Q  equals  5  in  a  fully  turbulent  region  [17]. 

The  correspondence  between  the  statistics  of  the  Monte  Carlo  particles  and  those  of  the  un¬ 
derlying  turbulent  reactive  flow  needs  careful  consideration  [1, 3].  A  primary  aim  of  the  modeling 
is  for  the  one-point,  one-time  joint  PDF  of  the  particles  to  accurately  represent  the  same  joint 
PDF  of  the  fluid  in  the  reactive  flow.  On  the  other  hand,  two-point,  one-time  statistics  are  radically 
different:  in  the  particle  system  the  properties  at  two  points  are  statistically  independent  (in  the 
infinite  particle  limit),  and  indeed  two  particles  may  have  the  same  location  x*  but  completely 
different  properties.  Of  particular  relevance  in  the  present  study,  is  the  question  of  correspon¬ 
dence  of  Lagrangian  statistics.  The  Langevin  equation  model  for  velocity  (Equation  4)  has  been 
constructed  to  be  consistent  with  the  Lagrangian  velocity  autocorrelation.  However,  the  mix¬ 
ing  models  have  been  developed  based  solely  on  one-time  Eulerian  statistics,  and  the  extent  to 
which  they  represent  Lagrangian  statistics  has  not  been  evaluated  even  in  simple  non-reactive 
flows. 

The  solution  procedure  for  the  above  equations  is  implemented  in  the  code  HYB2D  which 
implements  the  consistent  FV/particle  algorithm  [34].  The  HYB2D  code  has  been  fully  tested 
and  validated  in  various  papers,  e.g.  [23,  3 1, 34,  35].  This  work  slightly  extends  the  HYB2D  code 
to  output  the  Lagrangian  data. 

A  non-uniform  mesh  is  used  for  both  the  FV  solver  and  the  particle  tracking.  The  quantities 
at  the  mesh  level  are  interpolated  onto  particles  as  needed.The  statistics  of  particle  properties 
on  the  mesh  are  formed  from  particles  associated  with  the  cell.  The  mixing  between  particles  is 
performed  within  each  grid  cell. 

The  flames  we  are  interested  in  are  statistically  stationary.  The  solution  procedure  implemented 
in  HYB2D  is  a  pseudo-time  marching  procedure.  (A  local  time-stepping  algorithm  [43]  is  used 
for  the  marching  procedure.)  Starting  from  a  “reasonable”  initial  condition,  we  march  in  time 
steps  until  a  statistically  stationary  state  is  achieved.  We  continue  the  calculation  for  further  time 
steps  in  order  to  output  quantities  of  interest.  At  any  time  after  the  statistically  stationary  state 
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Figure  3.  The  computational  domain  and  the  tracking  domain  (consisting  of  four  sides  I,  II,  III,  and  IV) 
for  the  jet  flames. 


is  reached,  the  joint  PDF  is  represented  by  a  set  of  particles  which  (to  some  extent)  model  fluid 
particles. 

Conventionally,  in  the  PDF  calculations  only  Eulerian  data  are  output  for  postprocessing,  i.e., 
the  data  at  a  particular  time  step  after  the  statistically  stationary  state  has  been  reached.  Due  to 
the  Lagrangian  nature  of  the  numerical  method,  it  is  a  simple  matter  to  explore  the  Lagrangian 
data  by  simply  exporting  the  particle  data  for  many  time  steps  after  the  stationary  state  has  been 
reached,  so  that  the  Lagrangian  particle  trajectories  can  be  formed. 

Once  the  statistically  stationary  state  is  reached,  we  track  a  representative  number  Nt  of 
particles  through  the  active  part  of  the  flame  defined  as  the  tracking  domain  as  shown  in 
Figure  3.  The  geometry  of  the  computational  domain  and  the  tracking  domain,  and  the  val¬ 
ues  of  Nt  are  shown  in  Table  2  for  the  Sandia  flame  E  and  the  Cabra  lifted  H2/N2  jet  flame. 

In  the  particle  method,  there  is  a  particle  cloning  and  clustering  algorithm  designed  to  maintain 
an  approximately  uniform  number  of  particles  per  cell.  This  algorithm  creates  some  complications 
for  particle  tracking.  When  a  tracked  particle  is  cloned,  it  splits  into  two  or  more  (initially) 
identical  particles  of  less  weight.  We  arbitrarily  select  just  one  of  the  clones  to  continue  the 
particle  trajectory.  When  several  light  particles  are  clustered  to  form  one  heavier  particle,  the 
initial  identities  are  lost.  To  prevent  this  problem,  we  suppress  clustering  of  tracked  particles  (at 
a  small  cost  in  computational  accuracy  and  efficiency). 

Examination  of  the  particle  trajectories  in  physical  space  revealed  some  problems  with  the 
velocity-frequency  model  and  its  numerical  implementation.  These  are  discussed  in  the  Appendix 
where  a  method  of  alleviating  the  problem  is  described. 


Table  2.  Geometry  of  the  computational  domain  and  the  tracking  domain,  and  the  number  of  tracked 
particles  N,  for  the  Sandia  flame  E  and  the  Cabra  lifted  H2/N2  jet  flame.  See  Figure  3  for  definition  of 
locations. 


xu/D 

xA/D 

xa/D 

n/D 

rJD 

Nt 

Sandia  flame  E 

3.0 

45 

80 

10 

20 

2000 

Cabra  lifted  flame 

3.0 

30 

50 

10 

15 

2000 
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4.  Particle  trajectories  in  Sandia  flame  E 

PDF  calculations  of  the  Sandia  flame  E  are  performed  by  using  the  three  mixing  models,  EMST, 
IEM  and  modified  Curl.  The  Lagrangian  tracking  of  particles  is  conducted  to  investigate  the 
roles  of  reaction  and  mixing  in  the  regions  of  local  extinction  and  re-ignition.  We  focus  on  the 
particle  behavior  based  on  the  evolution  of  the  particle  temperature  in  the  mixture  fraction  space. 
As  in  the  scatter  plots  of  particles  in  Figures  1  and  2,  and  similar  to  the  DNS  analysis  in  [8], 
we  divide  the  particle  trajectories  into  two  groups:  continuous  burning  and  local  extinction.  For 
continuous  burning,  the  whole  particle  trajectory  remains  within  the  burning  region  (i.e.,  above 
the  “extinction  line”).  For  local  extinction,  some  segment  of  the  particle  trajectory  lies  in  the 
extinction  region  (below  the  “extinction  line”).  Physically,  continuous  burning  corresponds  to 
a  stretched  and  distorted  yet  still  continuous  non-premixed  laminar  flame  front,  and  the  local 
extinction  produces  holes  in  the  flame  front  [44,  45].  It  should  be  appreciated,  however,  that  in 
PDF  methods  there  is  no  representation  of  the  instantaneous  flame  structure.  For  ease  of  analysis, 
in  each  group,  we  further  sub-divide  the  particles  into  different  categories  based  upon  their 
mixture  fraction  at  the  trajectory’s  initial  position  xu  (see  Figure  3  and  Table  2  for  details),  i.e., 
fuel  region  (£  <  0.1),  oxidizer  region  (£  >  0.9),  pilot  stream  region  (0.22  <  f  <  0.55),  and  the 
intermediate  region  (all  other  values  of  f ). 

In  the  flames  considered  here,  the  particle  axial  distance  **(f)  is  an  increasing  function  of 
time  r.  To  some  extent,  the  axial  distance  can  be  viewed  as  a  time-like  variable  since  the  particles 
do  not  flow  backwards  in  the  axial  direction.  Since  we  are  more  interested  in  local  extinction  and 
re-ignition  at  different  axial  locations,  it  is  more  revealing  to  explore  the  particle  time  series  with 
respect  to  the  axial  distance  x/D  rather  than  with  respect  to  time. 


4. 1.  Trajectories  of  continuously  burning  particles 

The  trajectories  of  continuously  burning  particles  from  the  fuel  region  in  flame  E  are  shown  in 
Figures  4-6  for  PDF  calculations  using  the  different  mixing  models.  Only  25  particles  randomly 
chosen  from  the  tracking  dataset  are  shown  for  each  category.  The  circles  in  the  plots  show 
the  current  compositions  of  particles,  and  the  lines  connect  their  past  compositions.  For  each 
figure  of  particle  trajectories  (such  as  Figures  4—6),  the  supplementary  material  includes  a  cor¬ 
responding  animation.  These  animations  show  the  evolution  with  axial  distance  of  all  tracked 
particles’  compositions.  In  the  temperature-mixture  fraction  2-D  plane  in  Figures  4—6,  chemistry 
can  only  change  the  particle  positions  vertically  due  to  element  conservation  during  reaction 
(conservation  of  mixture  fraction  £),  while  mixing  can  move  the  particles  both  vertically  and 
horizontally. 

The  general  observations  on  the  trajectories  of  the  continuously  burning  particles  calculated 
using  the  different  mixing  models  are  the  following.  First,  the  mixture  fraction  of  particles  can 
vary  in  the  whole  mixture  fraction  range,  e.g.,  initially  is  greater  than  0.9  for  all  particles, 
while  later  £*  is  less  than  0. 1  for  some  particles  in  Figures  4-6.  As  discussed  before,  this  change 
is  solely  caused  by  mixing,  indicating  the  important  role  of  mixing  in  turbulent  combustion. 
Second,  different  particles  have  completely  different  trajectories  as  expected  in  a  turbulent  flow. 
Third,  at  different  stages  of  the  particle  evolution,  the  roles  of  reaction  and  mixing  are  different.  In 
Figures  4-6,  the  particles  from  the  fuel  region  tend  to  come  close  to  the  extinction  line  when  first 
approaching  the  stoichiometric  condition  (f  =  0.351).  At  around  the  stoichiometric  condition  in 
Figures  4-6,  we  can  observe  that  some  particles  suddenly  shoot  upward  (e.g.,  at  */D=18,  some 
particle  trajectories  become  nearly  vertical).  Apparently,  in  this  stage,  the  reaction  time  scale  is 
much  less  than  the  mixing  time  scale,  and  reaction  becomes  dominant.  A  quantitative  presentation 
of  the  relative  roles  of  mixing  and  reaction  is  discussed  in  Section  6. 
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Figure  4.  The  continuous  burning  particle  trajectories  from  the  fuel  region  in  flame  E  by  the  EMST  model. 
(This  [link]  provides  an  animation  of  these  particle  trajectories.) 
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Figure  5.  The  continuous  burning  particle  trajectories  from  the  fuel  region  in  flame  E  by  the  I  EM  model. 
(This  [link]  provides  an  animation  of  these  particle  trajectories.) 
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Figure  6  The  continuous  burning  particle  trajectories  from  the  fuel  region  in  flame  E  by  the  modified  Curl 
model.  (This  [link]  prov  ides  an  animation  of  these  particle  trajectories.) 


The  particle  behavior  simulated  by  the  different  mixing  models  is  qualitatively  different.  The 
particle  trajectories  produced  by  the  EMST  model  are  continuous  but  non-differentiable  [21]. 
These  trajectories  are  not  smooth  in  Figure  4.  The  trajectories  by  the  IEM  model  are  continuous 
and  differentiable,  and  the  simulated  trajectories  are  smooth  and  are  clear  in  Figure  5.  It  is 
easy  to  follow  each  particle  from  the  plot,  making  the  observation  and  analysis  much  easier. 
The  trajectories  arising  from  the  modified  Curl  model  are  discontinuous.  The  particles  jump 
in  composition  space,  possibly  resulting  in  the  direct  mixing  of  a  cold  fuel  particle  and  a  cold 
oxidizer  particle.  This  can  be  observed  at  x/D  =  45  in  Figure  6.  Two  particles  with  very  lean  and 
very  rich  mixtures  are  connected,  indicating  the  jumping  of  the  particle  from  the  one  side  to  the 
other  side  instantaneously.  (For  continuously  burning  particles,  by  definition,  their  compositions 
at  no  time  lie  in  the  extinction  region,  which  means  the  straight  lines  across  the  extinction  region 
at  x/D  —  45  in  Figure  6  correspond  to  an  instantaneous  jump  in  particle  composition.)  This  jump 
behavior  by  the  modified  Curl  makes  the  particle  trajectories  difficult  to  follow. 

Similar  observations  can  be  made  from  the  trajectories  of  continuously  burning  particles  from 
other  categories  (oxidizer  region,  pilot  stream  region  and  the  intermediate  region).  Due  to  space 
limitations,  these  trajectories  are  not  shown  here. 


4.2.  Trajectories  of  locally  extinguished  particles 

4.2  1  Particle  trajectories  from  the  EMST  model 

The  locally  extinguished  particle  trajectories  originating  from  different  regions  in  calculations 
using  the  EMST  mixing  model  are  shown  in  Figures  7-9. 
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Figure  7.  Trajectories  of  loeally  extinguished  partieles  from  the  fuel  region  in  flame  E  by  the  EMST  model. 
(This  [link]  provides  an  animation  of  these  particle  trajectories.) 
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Figure  8.  Trajectories  of  locally  extinguished  particles  from  the  pilot  stream  region  in  flame  E  by  the 
EMST  model.  (This  [link]  provides  an  animation  of  these  particle  trajectories.) 
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Figure  9.  Trajectories  of  locally  extinguished  particles  from  the  intermediate  region  in  flame  E  by  the 
EMST  model  .  (This  [link]  prov  ides  an  animation  of  these  particle  trajectories.) 


Figure  7  shows  the  trajectories  of  particles  initially  from  the  fuel  region.  In  the  temperature- 
mixture  fraction  space,  the  mixture  fraction  of  particles  first  decreases,  and  the  particles  come 
close  to  the  extinction  line.  Around  £  =  0.6,  the  particles  enter  the  extinction  region,  and  become 
locally  extinguished  according  to  our  criterion  for  extinction.  The  particle  trajectories  inside 
the  extinction  region  become  nearly  horizontal  (little  temperature  rise),  implying  that  mixing  is 
at  least  as  rapid  as  reaction.  The  return  of  the  particles  to  the  burning  region  corresponds  to 
re-ignition.  From  Figure  7,  two  different  re-ignition  processes  can  be  observed.  First,  at  around 
the  stoichiometric  condition,  some  trajectories  of  extinguished  particles  turn  and  move  upward 
to  return  to  the  burning  state,  e.g.,  at  x/D  =  18  in  Figure  7  we  can  observe  four  trajectories 
of  extinguished  particles  moving  dominantly  upward  to  return  to  the  burning  state.  During  this 
re-ignition  process,  the  mixture  fraction  of  the  particles  changes  slightly  while  the  temperature 
rises  by  more  than  600  K.  Although  re-ignition  is  the  result  of  mixing  and  reaction,  reaction 
seems  dominant  in  this  re-ignition  process,  which  is  similar  to  auto-ignition.  The  local  extinction 
induced  by  the  mixing  causes  the  coexistence  of  fuel  and  oxidizer  in  the  same  particle,  and  the 
temperature  of  these  particles  is  greater  than  1000  K,  e.g.,  one  particle  in  the  extinction  region  at 
x/D  =  15  in  Figure  7.  Given  appropriate  conditions  (e.g.,  induction  period  and  a  relatively  long 
mixing  time  scale),  the  auto-ignition  brings  the  particles  back  to  the  burning  state.  We  refer  to  this 
re-ignition  mechanism  as  an  auto-ignition  mechanism.  Second,  instead  of  auto-ignition,  the  other 
extinguished  particles  keep  moving  in  the  same  direction  (nearly  horizontally  to  the  left)  and 
re-enter  the  burning  state  on  the  lean  side  of  stoichiometric,  e.g.,  from  x/D  =  20  to  45.  During 
this  process,  mixing  is  at  least  as  rapid  as  reaction.  We  call  this  mechanism  the  mixing-reaction 
mechanism.  It  is  worth  mentioning  that  the  above  two  re-ignition  mechanisms  are  identified  in 
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the  DNS  study  [11].  The  particles  from  the  oxidizer  region  (not  shown)  behave  similarly  to  the 
particles  from  the  fuel  region.  The  two  different  re-ignition  processes  are  also  observed  there. 

Figure  8  shows  the  trajectories  of  the  particles  from  the  pilot  stream  region  in  flame  E  when 
using  the  EM  ST  model.  The  pilot  stream  is  used  to  stabilize  the  flame.  From  Figure  8,  for  this 
subset  of  particles  (all  of  which  enter  the  extinction  region  at  some  time),  they  dominantly  enter 
the  extinction  region  from  the  lean  and  rich  sides,  after  there  has  been  mixing  essentially  along 
the  fully  burnt  line,  Only  one  or  two  particles  enter  the  extinction  region  from  above  around 
stoichiometric. 

In  Figure  9  are  shown  the  trajectories  of  the  particles  from  the  intermediate  region  in  flame  E 
using  the  EMST  model.  Initially  (x/D  ^  7.5)  the  particle  composition  changes  due  to  mixing  and 
reaction  and  remains,  predominantly,  in  the  burning  region.  The  local  extinction  observable  for 
12  ^  x/D  ^  25  occurs  dominantly  by  mixing  drawing  rich  and  lean  particles  nearly  horizontally 
into  the  extinction  region. 

In  summary,  the  local  extinction  and  re-ignition  processes  in  the  Sandia  flame  E  are  illustrated 
by  tracking  particles  using  the  EMST  model.  Two  different  re-ignition  mechanisms  are  observed 
in  the  flame  by  using  the  EMST  model,  i.e.  auto- ignition  and  mixing-reaction.  The  investigation  of 
the  mixing  models  in  conjunction  with  the  large  eddy  simulations  using  the  DNS  data  by  Mitarai 
et  al.  [10]  demonstrates  the  very  good  performance  of  the  EMST  mixing  model  in  predicting  the 
particle  behavior  in  the  regions  of  local  extinction  and  re-ignition.  The  above  observed  particle 
behavior  by  the  EMST  is  expected  to  represent  the  actual  situation  qualitatively. 


4.2.2.  Particle  trajectories  using  the  IEM  and  modified  Curl  models 

The  trajectories  of  the  particles  from  the  fuel  region  in  flame  E  using  the  IEM  model  are  shown 
in  Figure  10.  Following  each  particle  trajectory,  we  can  observe  the  similar  local  extinction  and 
re- ignition  processes  as  in  the  case  of  the  EMST  model  (shown  in  Figure  7).  The  two  re-ignition 
mechanisms  can  also  be  identified:  auto-ignition  and  mixing-reaction.  The  re-ignition  process  for 
the  IEM  model,  however,  is  somewhat  different  from  that  for  the  EMST  model  in  Figure  7.  In 
Figure  7,  the  re-igniting  particles  for  the  EMST  model  tend  to  move  almost  horizontally  first  with 
a  slight  temperature  nse,  and  then  either  move  upward  due  to  the  auto-ignition  mechanism  or 
keep  moving  horizontally  due  to  the  mixing-reaction  mechanism  without  an  obvious  temperature 
drop  before  entenng  the  burning  region.  In  Figure  10,  however,  some  of  the  re-igniting  particles 
induced  by  the  auto-ignition  mechanism  experience  a  temperature  drop  before  ignition.  Almost 
all  the  re-igniting  particles  induced  by  the  mixing-reaction  mechanism  tend  to  decrease  their 
temperature  significantly  before  entering  the  burning  region.  For  the  IEM  model  at  the  early 
stages  of  re-ignition,  this  behavior  is  explained  by  the  fact  that  the  mean  temperature  (to  which 
the  particle  temperature  relaxes)  is  lower  than  that  of  the  particles,  which  are  about  to  re-ignite, 
since  these  are  the  hottest  particles  in  the  ensemble.  To  some  extent  this  reflects  the  physics  of 
the  problem  in  that  conduction  cools  fluid  at  a  local  temperature  maximum.  The  trajectories  of 
the  particles  from  the  other  regions  in  flame  E  using  the  IEM  model  (not  shown)  show  similar 
behavior  to  those  from  the  fuel  region  in  Figure  10. 

Figure  1 1  shows  the  particle  trajectories  from  the  fuel  region  using  the  modified  Curl  model. 
The  jumps  in  the  particle  properties  make  the  understanding  of  particle  behavior  more  difficult. 
The  particle  evolution  is  generally  quite  similar  to  the  IEM  model.  Local  extinction  and  re¬ 
ignition  are  predicted,  and  the  two  re-ignition  mechanisms  can  be  observed.  However,  similar  to 
the  IEM  model,  the  re-igniting  particles  tend  to  have  some  temperature  drop  before  or  during  the 
re-ignition,  which  is  not  observed  in  the  EMST  results.  The  trajectories  of  the  particles  from  the 
other  regions  in  flame  E  using  the  modified  Curl  model  (not  shown)  show  similar  behavior  to 
those  from  the  fuel  region  in  Figure  1 1 . 
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Figure  10.  Trajectories  of  locally  extinguished  particles  from  the  fuel  region  in  flame  E  using  the  1EM 
model.  (This  [link]  provides  an  animation  of  these  particle  trajectories.) 


Figure  1 1 .  Trajectories  of  locally  extinguished  particles  from  the  fuel  region  in  flame  E  using  the  modified 
Curl  model.  (This  [link]  provides  an  animation  of  these  particle  trajectories.) 
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In  this  sub-section,  the  local  extinction  and  re-ignition  processes  in  the  Sandia  flame  E  are 
illustrated  by  tracking  particles  using  the  IEM  and  modified  Curl  models.  The  two  re-ignition 
mechanisms  (auto-ignition  and  mixing-reaction)  identified  by  using  the  EMST  model  in  Figure 
7  are  also  observed  here  by  using  these  two  mixing  models.  However,  the  re-igniting  particles  by 
these  two  mixing  models  have  somewhat  different  behavior  from  those  by  the  EMST  model,  i.e., 
a  temperature  drop  before  or  during  the  re-ignition.  This  difference  in  re-ignition  by  the  different 
mixing  models  is  not  clear  yet  because  there  is  no  experimental  data  on  Lagrangian  trajectories 
in  the  flame.  Nevertheless,  the  two  identified  re-ignition  mechanisms  and  the  different  particle 
behavior  during  the  re-ignition  by  the  different  mixing  models  cannot  be  observed  with  the 
Eulerian  particle  data,  and  the  observations  contribute  to  our  understanding  of  the  performance 
of  the  models. 


5.  Particle  trajectories  in  Cabra  H2/N2  lifted  flame 

Previous  studies  [31,  32]  suggest  that  auto-ignition  is  a  dominant  mechanism  in  the  stabilization 
of  the  Cabra  H2/N2  lifted  flame.  In  addition,  both  experimentally  [33]  and  in  modeling  studies 
[31],  it  is  found  that  the  flames  are  extremely  sensitive  to  the  temperature  of  the  vitiated  coflow. 
To  further  understand  and  characterize  these  processes,  we  first  perform  auto-ignition  tests  in 
which  the  ignition  delay  time  (IDT)  is  calculated  as  a  function  of  mixture  fraction  and  coflow 
temperature.  The  initial  condition  of  the  tests  satisfies  Equations  (l}-( 2).  The  fuel  stream  and  the 
oxidizer  stream  in  the  tests  are  the  same  as  those  in  the  Cabra  lifted  flame.  The  coflow  (oxidizer) 
temperature  varies  from  Tc  =  1022  K  to  1080K.  Figure  12  shows  the  IDTs  of  the  mixture  for 
different  coflow  temperatures.  The  IDT  is  defined  here  as  the  time  when  the  mixture  temperature 
reaches  the  mid-point  between  the  initial  temperature  and  the  equilibrium  temperature.  The  strong 
sensitivity  of  the  IDTs  to  the  coflow  temperature  is  evident  from  the  plot,  which  is  consistent  with 
the  findings  in  [31].  The  shortest  IDTs  occur  at  the  very  fuel-lean  region,  around  £  =  0.04.  The 
stoichiometric  condition  is  £  =  0.47  in  the  flame.  The  IDTs  varies  by  three  orders  of  magnitude 
over  the  range  of  mixture  fraction  shown. 
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Figure  12.  The  ignition  delay  time  (IDT)  of  H2/N2/O2  mixture  for  different  coflow  temperature  Tc. 
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The  PDF  calculations  of  the  Cabra  H2/N2  lifted  flame  are  performed  by  using  the  three  mixing 
models,  EMST,  IEM  and  modified  Curl.  The  Lagrangian  tracking  of  particles  is  conducted  to 
investigate  the  roles  of  reaction  and  mixing  in  the  flame.  The  tracking  details  are  shown  in  Table  2. 
As  in  the  analysis  of  flame  E  in  Section  4,  we  focus  on  the  particle  behavior  based  on  the  evolution 
of  the  particle  temperature  in  mixture  fraction  space.  The  particle  trajectories  are  divided  into 
different  categones  based  upon  their  mixture  fraction  at  the  trajectory  initial  position  jcu,  i.e.,  fuel 
region  (£  <  0.1),  oxidizer  region  (£  >  0.9),  and  the  intermediate  region  between  the  fuel  and 
the  oxidizer  region.  For  each  category,  100  particles  randomly  chosen  from  the  tracking  dataset 
are  shown  in  the  following  figures,  whereas  all  tracked  particles  are  shown  in  the  corresponding 
animations  in  the  supplementary  material. 

Figure  13  shows  the  trajectories  of  the  particles  from  the  fuel  region  in  the  Cabra  lifted 
flame  using  the  EMST  model.  Initially,  (x/D  ^  9)  the  particles  move  in  the  plane  exclusively  by 
mixing.  A  particle  trajectory'  due  to  pure  mixing  is  a  nearly  straight  line  between  the  cold  fuel 
temperature  and  the  hot  coflow  temperature.  Pure  mixing  yields  a  partially  premixed  mixture 
of  fuel  and  oxidizer  at  different  mixture  fractions.  At  about  x/D  —  10,  some  particles  near  the 
oxidizer  side  start  to  ignite  first  due  to  their  short  IDTs  as  shown  in  Figure  12.  The  ignition 
mechanism  of  the  first  few'  particles  is  expected  to  be  auto-ignition,  similar  to  the  auto-ignition 
of  the  homogeneous  mixture  in  Figure  12. 

After  the  rapid  auto-ignition  of  the  first  few  particles,  these  relatively  hot  burnt  particles 
at  x/D  >  11  in  Figure  13  mix  with  adjacent  particles  in  composition  space,  thus  raising  their 
temperature  (and  radical  concentration)  and  hence  promoting  their  auto-ignition.  Therefore  the 
ignition  progressively  moves  to  richer  mixtures.  This  burning  process  is  not  exclusively  the 
auto-ignition  of  the  particles.  Both  reaction  and  mixing  play  important  roles.  A  plausible  physical 
picture  of  the  processes  involved  in  the  Cabra  flame  is  that  some  regions  under  the  fuel-lean 
condition  ignite  first  after  the  induction  period  given  an  appropnate  mixing  condition.  These 
ignition  spots  are  distributed  in  the  physical  space  separately.  After  the  high  temperature  ignition 
spots  are  formed,  they  propagate  toward  each  other  and  merge  into  a  connected  premixed  flame 
front.  This  picture  is  supported  by  the  DNS  study  of  the  auto-ignition  of  mixing  layers  between 
cold  fuel  and  hot  oxidizer  in  an  isotropic  and  homogeneous  turbulence  flow  [46].  The  evolution 
of  the  particles  using  EMST  in  Figure  13  is  consistent  with  this  picture,  even  though  the  spatial 
structure  of  the  instantaneous  flame  is  not  explicitly  represented.  We  simply  name  this  ignition 
process  as  mixing-ignition.  By  x/D  =  30  in  Figure  13,  all  the  particles  shown  reach  the  full  burnt 
state  close  to  the  equilibrium  line.  From  the  particle  trajectories  in  the  Cabra  lifted  jet  flame  using 
the  EMST  model,  we  can  observe  the  whole  mixing-reaction  process.  Four  stages  of  combustion 
can  be  identified,  i.e.  pure  mixing,  auto-ignition,  mixing-ignition,  and  fully  burnt.  Apparently, 
this  combustion  detail  cannot  be  observed  from  the  Eulerian  data  like  the  scatter  plot  in  Figure  2. 
The  particles  from  the  other  regions  in  the  Cabra  lifted  flame  by  using  EMST  model  (not  shown) 
show  the  similar  ignition  dynamics. 

In  Figure  14  are  shown  the  trajectories  of  the  particles  from  the  fuel  region  in  the  Cabra 
lifted  flame  using  the  IEM  model.  Pure  mixing  occurs  for  x/D  <  10.  At  the  locations  between 
x/D  =  10  and  1 1,  a  few  particles  near  the  oxidizer  side  (brought  there  by  mixing)  start  to 
auto-ignite.  After  the  auto-ignition  of  the  first  fewr  particles,  the  temperature  of  other  particles  in 
the  rich  region  is  raised  through  their  mixing  with  the  elevated  mean,  and  the  ignition  of  these 
particles  is  promoted.  Mixing  and  reaction  play  important  roles  in  this  ignition  process.  As  in 
the  EMST  model,  this  ignition  process  can  be  named  as  mixing-ignition.  The  particle  behavior 
in  this  ignition  process  for  the  IEM  model  is  slightly  different  from  that  for  the  EMST  model  in 
Figure  13.  The  EMST  model  is  local  in  composition  space.  Hence  the  burnt  particles  at  given  £* 
mix  with  particles  around  the  same  value  of  £*.  While  in  the  IEM  model  all  particles  mix  towards 
the  mean.  In  spite  of  the  different  particle  behavior,  four  stages  of  combustion  can  be  identified 
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Figure  13.  Particle  trajectories  from  the  fuel  region  in  the  Cabra  lifted  flame  using  the  EMST  model.  (This 
[link]  provides  an  animation  of  these  particle  trajectories.) 


Figure  14.  Particle  trajectories  from  the  fuel  region  in  the  Cabra  lifted  flame  using  the  IEM  model.  (This 
[link]  provides  an  animation  of  these  particle  trajectories.) 
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Figure  15.  The  particle  trajectories  from  the  fuel  region  in  the  Cabra  lifted  flame  by  the  modified  Curl 
model.  (This  [link]  prov  ides  an  animation  of  these  particle  trajectories.) 


for  the  IEM  as  those  in  the  EMST  model.  The  trajectories  of  the  particles  from  the  other  regions 
(not  shown)  show  similar  particle  behavior  for  IEM  model  as  in  Figure  14. 

In  Figure  15  are  shown  the  particle  trajectories  from  the  fuel  region  in  the  Cabra  lifted  flame 
using  the  modified  Curl  model.  The  modified  Curl  model  can  reproduce  the  same  four  combustion 
stages  in  the  flame  identified  by  the  previous  two  mixing  models.  As  in  the  IEM  model,  the  particle 
behavior  in  the  mixing-ignition  stage  by  the  modified  Curl  model  is  also  different  from  that  by 
the  EMST  model  due  to  the  non-localness  of  the  model  in  the  composition  space.  The  particle 
trajectories  from  the  other  regions  in  the  Cabra  lifted  flame  using  the  modified  Curl  model  (not 
shown)  show  the  same  behavior. 

In  this  section,  the  particle  trajectories  using  the  different  mixing  models  are  investigated 
in  the  Cabra  H2/N2  lifted  flame.  The  particle  behavior  by  the  IEM  and  modified  Curl  mod¬ 
els  is  different  from  that  by  the  EMST  model,  because  of  the  non-localness  of  the  IEM  and 
modified  Curl  models  compared  to  the  localness  property  of  the  EMST  model.  In  spite  of  the 
different  individual  particle  behavior,  the  overall  combustion  processes  revealed  by  the  differ¬ 
ent  mixing  models  are  similar,  and  four  stages  of  combustion  in  the  flame  can  be  identified, 
i.e.  pure  mixing,  auto-ignition,  mixing-ignition,  and  fully  burnt.  In  some  sense,  this  finding  is 
consistent  with  the  DNS  study  of  an  auto-ignition  problem  in  homogeneous  isotropic  turbu¬ 
lence  [46],  even  though  the  particles  do  not  provide  a  direct  representation  of  spatial  struc¬ 
ture.  This  contributes  to  our  understanding  of  the  model  performance  in  the  turbulent  lifted  jet 
flames. 
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Figure  16.  Trajectories  (up  to  x/l)  =  45)  color-coded  by  parameter  r\  (Equation  10)  of  locally  extinguished 
particles  from  the  fuel  region  in  flame  E  using  the  IEM  model. 


6.  Roles  of  mixing  and  reaction  during  re-ignition  and  auto-ignition 

In  the  previous  sections,  the  roles  of  mixing  and  reaction  during  the  particle  evolution  are  discussed 
qualitatively.  The  relative  importance  of  mixing  and  reaction  for  each  particle  as  it  evolves  can 
be  quantified  by  examining  the  mixing  rate  given  by  the  mixing  models  relative  to  the  reaction 
rate.  For  this  purpose  we  define  the  mixing  rate  M  and  reaction  rate  S  for  each  particle  as, 


mix 


Figure  1 7.  Particle  trajectories  (up  to  x/ D  =  30)  color-coded  by  parameter  rj  (Equation  10)  from  the  fuel 
region  in  the  Cabra  lifted  flame  using  the  IEM  model 
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and 


S  = 


1  dT 
Tref  dt 


(9) 


where  Tref  —  2000  K,  and  in  Equations  (8)  and  (9)  the  rates  of  change  pertain  solely  to  the  effects 
of  mixing  and  reaction,  respectively.  (In  the  computations,  these  quantities  are  readily  evaluated 
based  on  the  particle  properties  before  and  after  the  mixing  and  reaction  fractional  steps.) 

The  relative  importance  of  mixing  and  reaction  can  be  quantified  by  the  parameter 


n  = 


5 

S  + M’ 


(10) 


which  varies  between  zero  (corresponding  to  no  reaction)  and  one  (corresponding  to  no  mixing). 

We  focus  our  investigation  on  the  IEM  mixing  model.  For  the  modified  Curl  model,  the 
composition  changes  discontinuously,  and  so  /dt  |mix  is  not  well  defined.  For  the  EMST  mixing 
model,  d£ /dt\m[x  exhibits  large  fluctuations  of  small  time  scale  which  obscure  the  picture. 

Figure  1 6  shows,  for  flame  E,  trajectories  color-coded  by  the  parameter  r]  of  locally  extin¬ 
guished  particles  initially  from  the  fuel  region.  From  the  figure,  we  can  clearly  see  the  relative 
importance  of  mixing  and  reaction  during  re-ignition.  For  the  particles  re-igniting  due  to  auto¬ 
ignition  mechanism,  the  reaction  is  dominant  (rf  ^  1.0)  when  the  particles  shoot  upward  at  around 
stoichiometric  condition.  For  those  particles  which  re-ignite  due  to  mixing-reaction,  either  both 
mixing  and  reaction  are  important  (e.g.,  rf  &  0.3  for  one  particle  entering  the  burning  region  at 
about  T  =  1 100K),  or  mixing  is  dominant  due  to  the  low  temperature. 

In  Figure  17,  the  particle  trajectories  originating  from  fuel  stream  in  the  Cabra  flame  are 
shown.  In  the  current  simulation  of  the  Cabra  flame,  the  initial  ignition  process  occurs  when 
particles  leave  the  pure  mixing  line  between  cold  fuel  and  hot  coflow.  From  the  figure,  it  may 
be  seen  that  the  particle  starting  to  ignite  at  fuel-lean  side  leave  the  pure  mixing  line  dominantly 
by  reaction,  corresponding  to  the  auto-ignition  identified  before.  The  particles  leaving  the  pure 
mixing  line  on  the  fuel-rich  side  expenence  two  stages:  a  mixing-dominant  stage  to  raise  the 
particle  temperature  to  about  1000K,  and  reaction-dominant  stage  to  raise  the  particle  temperature 
close  to  the  equilibrium.  Both  mixing  and  reaction  are  important  for  the  ignition  of  these  particles, 
and  they  correspond  to  the  previously  identified  mixing-ignition  process.  The  particle  trajectories 
along  w  ith  the  rates  of  mixing  and  reaction  provide  insights  on  the  roles  of  mixing  and  reaction 
during  re-ignition  and  auto-ignition. 


7.  Conclusion 

Lagrangian  PDF  investigations  of  the  Sandia  piloted  flame  E  and  the  Cabra  H2/N2  lifted  flame 
are  performed  to  help  obtain  a  deeper  understanding  of  the  modeling  of  local  extinction,  re- 
lgnition  and  auto-ignition  in  these  flames.  Eulerian  scatter  plots  (shown  in  Figures  1-2)  of  the  two 
flames  from  the  PDF  calculations  are  reviewed  to  show  the  limitations  of  one-time  statistics.  A 
Lagrangian  particle  tracking  procedure  is  implemented  in  the  code  HYB2D.  Lagrangian  particle 
data  are  extracted  from  the  PDF  calculations  after  the  statistically  stationary  state  is  reached,  in 
order  to  explore  the  PDF  result  more  comprehensively. 

Lagrangian  particle  tracking  in  the  PDF  calculations  of  Sandia  flame  E  is  performed  for  the 
different  mixing  models,  EMST,  IEM  and  modified  Curl.  The  particle  trajectories  are  divided  into 
two  groups,  continuous  burning  and  local  extinction.  For  each  group,  the  trajectories  are  further 
sub-divided  into  different  categories  based  on  the  original  particle  locations:  the  fuel  steam,  the 
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oxidizer  stream,  the  pilot  stream,  and  the  intermediate  region.  The  particle  trajectories  given  by  the 
different  mixing  models  are  different,  i.e.,  continuous  but  non-differentiable  by  EMST,  continuous 
and  differentiable  by  IEM,  and  discontinuous  by  modified  Curl.  All  three  mixing  models  reproduce 
the  local  extinction  and  re-ignition  processes  reasonably.  Two  different  re- ignition  mechanisms 
are  identified,  the  auto-ignition  mechanism  and  the  mixing-reaction  mechanism. 

Homogeneous  auto-ignition  tests  for  the  same  condition  as  in  the  Cabra  H2/N2  lifted  flame  are 
conducted.  The  lowest  ignition  delay  time  (IDT)  occurs  at  a  very  fuel-lean  condition  for  a  range 
of  coflow  (oxidizer)  temperatures.  The  strong  sensitivity  of  the  IDTs  to  the  coflow  temperature  is 
observed  which  is  also  reported  in  previous  PDF  calculations  of  the  Cabra  lifted  flame  [31]. 

Lagrangian  particle  tracking  in  the  PDF  calculations  of  the  Cabra  H2/N2  lifted  is  also  per¬ 
formed  for  the  different  mixing  models.  The  particle  trajectories  are  divided  into  different  cat¬ 
egories  based  on  the  original  particle  locations:  the  fuel  stream,  the  oxidizer  stream,  and  the 
intermediate  region.  The  models  reproduce  the  whole  auto-ignition  process  reasonably.  Four 
stages  of  combustion  in  the  Cabra  flame  are  identified  in  the  calculations,  i.e.,  pure  mixing, 
auto-ignition,  mixing-ignition,  and  fully  burnt. 

The  roles  of  mixing  and  reaction  during  re-ignition  and  auto-ignition  are  investigated  by  using 
IEM.  The  relative  importance  of  mixing  and  reaction  is  quantified  for  particles  during  re-ignition 
and  auto- ignition. 
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Appendix:  Particle  trajectories  in  physical  space 

Figure  18  shows  the  particle  trajectories  in  physical  space  for  the  Cabra  H2/N2  lifted  flame  and 
for  the  Sandia  flame  E,  from  calculations  using  the  EMST  mixing  model.  The  top  border  of 
the  plot  is  the  free-stream  boundary.  When  the  particles  from  the  turbulent  jet  approach  the 
non-turbulent  free  stream,  their  velocity  should  relax  rapidly  towards  the  free  stream  velocity. 
However,  from  Figure  18,  it  may  be  seen  that  some  particles  shoot  into  the  free  stream  with  little 
or  no  relaxation  of  velocity,  and  are  then  reflected  off  the  free-stream  boundary  (according  to 


Figure  18.  Particle  trajectories  in  physical  space  for  the  Cabra  H2/N2  lifted  jet  flame  and  for  the  Sandia 
flame  E  from  the  original  turbulence  frequency  model. 
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the  specified  boundary  condition).  This  non-physical  behavior  of  the  particles  is  not  expected 
to  occur  when  the  conditional  mean  frequency  Q  (Equation  (7))  is  used  to  define  the  time-scale 
in  the  stochastic  turbulence  frequency  model  (Equation  (5))  [17].  The  turbulence  frequency  in 
the  turbulent  region  is  much  greater  than  that  in  the  non-turbulent  free  stream.  Ideally,  if  one 
turbulent  fluid  particle  V  dives  into  one  mesh  cell  C  in  the  laminar  free  stream  environment,  the 
turbulence  frequency  of  particle  V  will  be  the  only  frequency  used  to  determine  £2  (=  Cq  *  co 
where  cop  is  the  turbulence  frequency  of  particle  V)  in  cell  C  according  to  Equation  (7).  From  the 
Langevin  model  Equation  (4),  the  velocity  of  particle  V  will  decay  rapidly  at  rate  Q  toward  the 
free  stream  velocity.  However,  under  certain  circumstances,  the  particle  V  is  not  the  only  particle 
chosen  to  determine  the  conditional  mean  frequency  Q.  In  the  current  situation  (of  a  cylindrical 
coordinate  system),  the  initial  mass  of  a  particle  is  linearly  proportional  to  the  particles  radial 
location.  Compared  to  the  particle  mass  in  the  free  stream,  the  mass  of  the  particle  V  originating 
closer  to  the  axis  is  much  less.  In  this  case,  the  mass- weighted  mean  frequency  5  in  cell  C  is  close 
to  the  mean  frequency  in  the  free  stream.  On  the  other  hand,  numerically  the  turbulence  frequency 
for  the  particles  in  the  free  stream  is  not  exactly  the  same,  i.e.,  there  are  small  fluctuations  in 
particle  turbulence  frequency  in  the  free  stream.  These  fluctuations  are  caused  by  the  Wiener 
process  in  Equation  (5)  and  by  the  disturbance  caused  by  the  turbulent  fluid  particle  V.  Assume 
that  the  maximum  of  the  turbulence  frequency  in  cell  C  is  co*m2  when  particle  V  is  excluded.  It 
may  happen  that  co*m 2  is  greater  than  5,  so  that  the  possibly  massive  particle  with  frequency 
is  included  in  the  calculation  of  Q.  Therefore  the  value  of  Q  is  much  less  than  the  desired  value 
Cq  •  (Op.  In  another  words,  the  value  of  Q  is  greatly  under-estimated,  and  so  is  the  decaying  rate 
of  the  velocity  of  the  particle  V. 

In  spite  of  the  non-physical  behavior  of  the  particles  described  above,  it  should  be  appreciated 
that  this  behavior  occurs  with  low  probability  (less  than  1%  of  the  number  of  tracked  particles 
in  the  Cabra  lifted  flame,  and  about  5%  in  the  Sandia  flame  E).  They  do  not  influence  the 
statistics  significantly.  Since  the  particle  behavior  is  relatively  more  important  in  this  work,  we 
try  to  eliminate  or  reduce  the  non-physical  behavior  in  the  following  ad  hoc  way.  In  the  above 
discussed  case,  the  particle  V  should  be  chosen  as  the  only  particle  to  determine  Q,  in  spite 
of  the  small  fluctuations  of  the  turbulence  frequency  in  the  free  stream.  If  we  can  identify  this 


x/D 

Figure  1 9.  Particle  trajectories  in  the  physical  space  for  the  Cabra  H2/N2  lifted  jet  flame  and  for  the  Sandia 
flame  E  with  an  ad  hoc  revision  to  the  turbulence  frequency  model. 
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particular  case,  we  can  avoid  this  problem  by  calculating  £2  using  particle  V  only.  We  use  the 
follow  ing  criteria  to  identify  the  case.  If  in  one  grid  cell,  the  maximum  turbulence  frequency  (Op 
of  a  particle  V  is  much  greater  than  the  Favre  mean  frequency  S  (o)p  >  c\  •  5),  and  also  much 
greater  than  the  frequency  of  all  other  particles  in  the  cell  (<Dp  >  ci  •  (o*m2 ),  and  the  mass  of  the 
particle  m p  is  much  less  than  the  average  mass  of  the  particles  (m)  {nip  <  C3  •  (m>),  we  then  use 
particle  V  exclusively  to  determine  £2.  The  constants  are  chosen  as  c\  =  15,  ci  =  8,  and  C3=0.2. 
The  particle  trajectories  in  physical  space  obtained  with  this  special  treatment  of  the  frequency 
model  are  shown  in  Figure  19.  The  non-physical  particles  disappear  in  the  Cabra  lifted  flame,  and 
the  number  of  the  non-physical  particles  is  significantly  reduced  in  Sandia  flame  E.  The  special 
treatment  improves  the  practical  performance  of  the  turbulence  frequency  model  to  some  extent. 
The  remaining  non-physical  particles  in  the  Sandia  flame  E  shown  in  Figure  19  are  due  to  the 
limitation  of  the  criteria.  Using  a  different  specification  of  Ci,  q,  and  C3,  we  can  eliminate  the 
non-physical  particles  completely,  but  it  also  then  affects  the  turbulence  region,  In  this  work, 
we  use  this  ad  hoc  revision  to  reduce  the  number  of  non-physical  particles.  Those  particles  not 
caught  by  the  criteria  are  removed  from  the  particle  subset  in  the  discussion.  This  ad  hoc  revision 
to  the  frequency  model  does  not  change  the  statistics  of  the  velocity  and  composition  fields.  We 
appreciate  that  the  criteria  does  not  guarantee  the  convergence  of  the  method  when  the  number 
of  particles  tends  to  infinity. 
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We  describe  a  particle  position  time  advancement  algorithm  that  is  designed  for  use  with 
several  subgrid  velocity  reconstruction  schemes  used  in  LES/FDF  methods,  and  potentially 
in  other  applications.  These  reconstruction  schemes  yield  a  subgrid  velocity  field  with 
desirable  divergence  properties,  but  also  with  discontinuities  across  cell  faces.  Therefore, 
a  conventional  time  advancement  algorithm,  such  as  second-order  Runge-Kutta  (RK2), 
does  not  perform  as  well  as  it  does  with  a  smooth  velocity  field.  The  algorithm  that  we 
describe,  called  Multi-Step  RK2  (MRK2),  builds  upon  RK2  by  breaking  up  the  time  step  into 
two  or  more  substeps  whenever  a  particle  crosses  one  or  more  velocity  discontinuities. 
When  used  in  conjunction  with  the  parabolic  edge  reconstruction  method.  MRK2  performs 
considerably  better  than  RK2:  both  the  final  position  of  an  advected  particle,  and  the  final 
area  of  a  2D  infinitesimal  area  element  are  second-order  accurate  in  time  (as  opposed  to 
first-order  accurate  for  RK2).  Furthermore,  MRK2  has  the  theoretical  advantage  that  it  bet¬ 
ter  preserves  the  continuity  of  the  mapping  between  initial  and  final  particle  positions. 

©  2008  Elsevier  Inc.  All  rights  reserved. 


1.  Introduction 

In  this  paper,  we  consider  some  aspects  of  particle  tracking  in  hybrid  finite-volume/particle  PDF  methods  for  turbulent 
reactive  flows.  In  these  methods,  and  potentially  in  other  CFD  applications,  we  have  a  large  number  of  particles  (on  the  order 
of  107)  whose  positions  are  initially  random  and  evolve  by  the  ODE: 

4£  =  ir  =  u(x*(r),o,  (i) 

where  X*(t)  is  a  particle’s  position  as  a  function  in  time,  and  U(X*(t),t)  is  the  velocity  experienced  by  the  particle.  Usually, 
this  velocity  consists  of  a  deterministic  component  and  a  random  term  which  is  part  of  the  turbulence  model  [5]:  in  the  pres¬ 
ent  paper,  we  consider  the  deterministic  part  only.  In  general,  velocity  information  is  available  at  discrete  locations  on  a  fi¬ 
nite-volume  (FV)  grid,  and  therefore  a  particle  tracking  algorithm  consists  of  two  parts:  a  velocity  interpolation  scheme 
which  interpolates  the  FV  velocity  onto  the  particle  locations,  and  a  time  advancement  scheme  which  updates  the  particle 
locations,  using  the  interpolated  velocity. 

In  recent  research  on  FV/particle  PDF  methods  for  reacting  flows,  accurate  particle  tracking  has  been  recognized  as  an 
important  condition  for  maintaining  numerical  consistency  between  the  finite-volume  and  the  particle  aspects  of  the  solu¬ 
tion.  Muradoglu  et  al.  [2]  define  mean  particle  mass  density,  q(x  t),  as  the  expectation  of  the  total  mass  of  particles  in  an 
infinitesimal  region,  normalized  by  that  region’s  volume: 

q  =  (m*<5(X*  -  x))  (2) 
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(where  rrv  denotes  a  particle’s  mass),  and  have  shown  that  an  important  consistency  condition  is  that  the  mean  particle  mass 
density  should  remain  equal  to  the  FV  mean  density: 

Q  =(P).  (3) 

provided  that  the  initial  conditions  are  consistent,  g(x.  0)  =  (p(x,  0)).  The  meaning  of  Eq.  (3)  is  that  the  expected  mass  of  all 
particles  inside  a  certain  region  S  must  be  equal  to  the  expected  mass  of  fluid  in  S  as  given  by  the  FV  density  field.  In  the 
incompressible  case,  for  example,  this  means  that  if  particles  are  initially  uniformly  distributed  in  the  computational  do¬ 
main,  then  they  should  remain  so  for  all  time. 

In  order  to  satisfy  the  above  consistency  condition,  Muradoglu  et  al.  [2]  employ  a  position-correction  algorithm  which 
introduces  a  small  displacement  in  the  particle  positions  after  each  time  step,  in  order  to  enforce  Eq.  (3).  A  similar  posi¬ 
tion-correction  algorithm  has  been  employed  by  Zhang  and  Haworth  [4]. 

In  a  different  approach,  Jenny  et  al.  [3]  note  that  if  particles  move  with  a  velocity: 

IT  =  U  +  IT,  (4) 

where  U  is  a  deterministic  component  interpolated  from  the  FV  velocity,  and  u*  is  a  random  component  with  zero  mean, 
then  the  following  equation  holds: 

(^■+  0  •  v)  'nq  =  -V  •  U,  (5) 

which  has  the  same  form  as  the  mean  continuity  equation: 

(!+U.v)ln(p)  =  -V.U,  (6) 

where  U  is  the  mean  velocity.  Therefore,  if  Eq.  (3)  is  satisfied  at  t  =  0,  it  will  also  be  satisfied  implicitly  at  future  time,  pro¬ 
vided  that  the  velocity  interpolation  scheme  yields  accurate  values  for  the  reconstructed  velocity  and  its  divergence,  and 
provided  that  an  accurate  time  advancement  scheme  is  used.  Addressing  the  velocity  interpolation  issue,  Jenny  et  al.  [3] 
introduce  a  2D  velocity  interpolation  scheme  with  desirable  divergence  properties.  McDermott  and  Pope  ( 1  ]  improve  upon 
this  scheme,  and  extend  it  to  3D,  calling  the  new  scheme  the  parabolic  edge  reconstruction  method  (PERM).  It  has  been 
shown  [  1  ]  that  PERM  performs  better  than  standard  multilinear  interpolation  in  satisfying  the  above-mentioned  consistency 
condition  in  the  particle  tracking  limit  (i.e.,  when  there  are  no  velocity  fluctuations). 

In  the  present  work,  we  make  a  further  improvement  in  particle  tracking  by  using  a  time-stepping  algorithm  which  has 
been  specifically  designed  for  use  in  conjunction  with  PERM.  The  new  algorithm,  called  Multi-Step  Runge-Kutta  2  (MRK2)  is 
quite  similar  to  the  standard  second-order  Runge-Kutta  (RK2),  but  it  provides  a  more  accurate  treatment  of  particles  which 
cross  a  velocity  discontinuity.  Although  MRK2  is  motivated  by  the  PERM  reconstruction,  it  can  also  be  used  as  an  alternative 
to  RK2  in  all  applications  in  which  particles  need  to  be  advected  though  a  velocity  field  with  discontinuities  at  known 
locations. 

To  illustrate  the  benefits  of  PERM  and  MRK2,  Fig.  1  shows  final  particle  distributions  for  a  simple  2D  incompressible  test 
flow  with  a  uniform  initial  particle  distribution.  In  this  case,  the  final  particle  distribution  should  be  uniform  as  well.  It  can  be 
seen  in  Fig.  1  that,  when  RK2  is  used  as  the  time  advancement  scheme,  the  PERM  velocity  interpolation  yields  a  more  uni¬ 
form  final  particle  distribution  than  standard  bilinear  interpolation.  We  can  also  see  that  the  final  particle  distribution  be¬ 
comes  even  more  uniform  when  the  new  time-stepping  scheme,  MRK2,  is  used  in  conjunction  with  the  PERM  velocity 
interpolation. 

1.1.  Properties  of  the  PERM  velocity  reconstruction 


The  reader  is  referred  to  [1]  for  a  thorough  description  of  PERM  and  its  properties  -  here  we  focus  only  on  its  aspects 
which  are  relevant  to  the  MRK2  time-stepping  scheme. 

The  PERM  velocity  reconstruction  scheme  uses  discretized  velocity  information  in  the  form  of  face-average  velocity  com¬ 
ponents  in  the  direction  normal  to  a  given  grid  cell  face  (we  shall  refer  to  this  as  the  ‘FV  velocity’,  for  the  sake  of  brevity).  The 
PERM  reconstructed  velocity  within  a  given  grid  cell  depends  only  on  the  FV  velocity  on  the  faces  of  that  cell  and  its  nearest 
neighboring  cells.  Its  functional  form,  in  two  dimensions  (the  extension  to  3D  is  rather  intuitive)  is 


u(q,r)  =  (1  -r) 


v(<?.r)  =  (1  -  q ) 


where  u,  v  are  respectively  the  horizontal  and  vertical  velocity  components,  N.  S,  E,  W  is  the  standard  compass  convention  for 
2D  grids,  q,  r  are  respectively  the  horizontal  and  vertical  local  coordinates  (e.g.,  q  =  0  on  the  West  face  and  q  =  1  on  the  East 
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Fig.  1.  Comparison  between  the  performance  of  the  bilinear  and  PERM  velocity  interpolation  schemes,  and  between  the  RK2  and  MRK2  time-stepping 
schemes.  We  consider  the  2D,  periodic  incompressible  flow  field  given  by  Eqs.  (16)  and  (17)  in  Section  3.  409,600  particles  are  initially  uniformly  spaced 
inside  the  domain  (top  left).  Each  particle  is  represented  by  a  light  grey  dot,  except  for  those  initially  in  the  bottom  left  corner,  which  are  dark  grey,  for 
visualization  purposes.  The  particles  are  advected  using  different  velocity  interpolation  and  time  integration  schemes.  Because  the  test  flow  is 
incompressible  and  the  particles  are  initially  uniformly  distributed,  they  should  remain  so.  Top  right,  final  particle  distribution  using  bilinear  interpolation 
on  a  4  x  4  grid  and  RK2  time  integration.  Bottom  left:  final  particle  distribution  using  PERM  interpolation  on  a  4  x  4  grid  and  RK2  time  integration.  Bottom 
right:  final  particle  distribution  using  PERM  interpolation  on  a  4  x  4  grid  and  MRK2  time  integration.  It  can  be  seen  that  the  combination  of  PERM 
interpolation  and  MRK2  time  stepping  performs  best  at  preserving  the  uniformity  of  the  particle  distribution. 


face.  Note  that  q  as  used  here  is  different  from  the  mean  particle  mass  density  defined  in  Eq.  (2)),  and  Us,  etc.  are 

parameters  which  are  determined  by  an  algorithm  described  in  |1). 

Some  of  the  important  properties  of  PERM  are  enumerated  below.  They  indicate  why  a  standard  time-stepping  scheme, 
such  as  RK2,  does  not  perform  well  with  this  particular  velocity  reconstruction.  For  simplicity,  we  consider  here  a  grid  com¬ 
posed  of  2D  or  3D  cubes,  of  side  Ax. 


1.  For  a  given  velocity  component,  in  the  component  direction,  the  reconstructed  field  is  continuous,  piecewise  parabolic, 
and  second-order  accurate  with  respect  to  Ax. 

2.  The  flux  of  the  reconstructed  velocity  field  through  any  given  grid  cell  face  is  identically  equal  to  the  flux  implied  by  the 
FV  face-normal  velocity. 

3.  For  a  given  velocity  component,  in  the  component-normal  directions,  the  reconstructed  field  is  piecewise  continuous, 
piecewise  linear,  and  second-order  accurate  with  respect  to  Ax. 
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4.  The  magnitudes  of  the  velocity  discontinuities  implied  by  Property  3  are  proportional  to  Ax2. 

5.  The  divergence  of  the  reconstructed  velocity  field  is  bi-  or  tri-hnear,  piecewise  continuous,  and  second-order  accurate 
with  respect  to  Ax. 

6.  The  magnitudes  of  the  divergence  discontinuities  implied  by  Property  5  are  proportional  to  Ax2. 

It  should  be  noted  that,  as  described  in  [1],  the  PERM  reconstruction  scheme  requires  a  Cartesian  grid  (either  uniform  or 
non-uniform  in  each  coordinate).  It  is  also  possible  to  extend  PERM  to  cylindrical  coordinate  grids. 

Furthermore,  it  should  be  noted  that  for  a  velocity  field  with  the  functional  form  of  PERM  it  is  not  possible  to  integrate  Eq. 
( 1 )  analytically  -  to  see  why  this  is  so,  note  that  any  bi-  or  tri-linear  velocity  field  is  a  particular  case  of  a  PERM  field,  and  Eq. 
(1 )  cannot  be  integrated  exactly  for  a  general  bi-  or  tri-linear  velocity  field,  even  if  it  is  time-independent  (the  differential 
equation  of  the  Lorentz  attractor  |6],  for  example,  is  a  particular  case  of  a  tri-linear  velocity  field).  Therefore,  we  have  to  uti¬ 
lize  a  numerical  time-integration  scheme  for  the  integration  of  Eq.  (1)  -  a  natural  candidate  is  second-order  Runge-Kutta 
(RIG). 

However,  second-order  time  accuracy  for  the  RK2  scheme  requires  that  the  velocity  field  is  everywhere  continuous,  and 
differentiable  (in  both  space  and  time)  everywhere  except  on  a  set  of  measure  zero.  Therefore,  due  to  the  discontinuities  in 
the  PERM  reconstructed  velocity,  an  RK2  solution  is  only  first-order  accurate  in  time,  for  a  fixed  grid  size.  Also,  applying  RK2 
to  a  discontinuous  velocity  field  leads  to  the  violation  of  an  important  continuity  principle,  namely:  if  two  particles  are  ini¬ 
tially  infinitesimally  close,  they  remain  so.  These  two  statements  are  demonstrated  in  the  next  section. 

On  the  other  hand,  Items  4  and  6  above  imply  that  for  a  fine  mesh  the  discontinuities  in  velocity  and  divergence,  and  their 
detrimental  effect  on  the  performance  of  RK2,  are  negligible.  More  specifically,  if  we  decrease  both  the  grid  size  and  time 
step,  keeping  the  Courant  number  fixed,  a  solution  which  combines  PERM  and  RK2  is  second-order  accurate.  Nevertheless, 
it  is  often  the  case  in  large  eddy  simulation/filtered  density  function  (LES/FDF)  methods  that  the  grid  used  provides  only  mar¬ 
ginal  resolution  of  the  filtered  velocity  field,  and  hence  the  discontinuities  in  the  PERM  reconstructed  velocity  field  should 
not  be  neglected. 


2.  Examination  of  RK2  for  discontinuous  velocity  fields 


Let  us  begin  with  a  simple  example  of  the  problems  that  we  encounter  when  we  use  RIG  on  a  discontinuous  velocity  field. 
Instead  of  thinking  about  a  grid,  consider  an  infinite,  time-independent  velocity  field  in  2D,  which  has  the  following  form: 


U(x,y)  =  1 

V(x,y)  =  \\ 


for  x  <  0 
for  x  ^  0 


(7) 

(8) 


Here,  we  can  think  of  the  regions  x  <  0  and  x  >  0  as  two  grid  cells,  with  the  region  x  0  as  the  face  between  them.  Note 
that  this  velocity  field  has  the  kind  of  discontinuity  which  can  be  caused  by  PERM:  U  is  continuous  in  the  x-direction,  and  V  is 
continuous  in  the  y-direction,  but  V  is  not  continuous  in  the  x-direction,  its  component-normal  direction.  Now,  consider  two 
particles  at  time  t  =  0:  one  at  (Xi.o,  Yi.o)  =  (-At  -  e,0)  and  the  other  at  (X20  Y2.o)  =  (-At  +  e,0),  where  0  <  e  «:  At.  An  RK2 
step  of  length  At  has  the  form: 

X(,)  =  Xo  +  AtU(Xo,0),  (9) 

X^2=Xo+iAt(U(Xo.0)  +  U,X"',At)).  (10) 

Substituting  their  initial  positions  for  Xo  in  Eqs.  (9)  and  (10),  we  determine  that  the  final  positions  of  the  two  particles, 
after  one  RK2  step  of  length  At,  are  (X*2i  Y*K42f)  =  (-e  0)  and  (X^,  Y^S)  =  (e,  At/2).  On  the  other  hand,  with  perfect  time 
advancement,  the  final  particle  positions  at  time  At  are  (Xi^,Yt.At)  =  (-e,0)  and  (X2Af,  Y2  ^)  =  (e,e/2).  Therefore,  we  can 
see  that  the  RIG  position  of  the  second  particle  differs  from  the  correct  position  by  a  distance  of  (At  -  e)/2.  Keeping  in  mind 
that  we  set  e  <  At,  we  have  that  this  single  time  step  introduces  an  O(At)  error  in  the  final  position  of  the  second  particle. 
Therefore,  RK2  is  first-order  accurate  in  time  overall,  and  for  a  time  step  that  crosses  a  discontinuity  it  is  zeroth-order  accu¬ 
rate  (in  the  sense  that  if  all  time  steps  introduced  an  error  of  the  same  magnitude  as  that  introduced  by  a  discontinuity- 
crossing  step,  the  overall  error  in  the  final  position  of  a  particle  would  be  0(Af°)). 

Furthermore,  if  we  take  the  limit  as  e  —  0,  we  can  see  that  the  distance  between  the  initial  positions  of  the  two  particles 
goes  to  zero,  but  the  distance  between  the  final  positions  after  the  RK2  step  goes  to  At/2.  As  previously  mentioned,  this  is  a 
violation  of  a  fundamental  continuity  principle,  namely  that  if  X(Y,  t)  is  the  Lagrangian  mapping  between  a  particle’s  initial 
and  final  positions  (i.e.,  if  a  particle  is  at  Y  at  time  0,  it  will  be  at  X(Y,  t)  at  time  t),  then  X(Y,  t)  is  continuous  in  both  Y  and  t, 
except  at  regions  where  there  is  a  velocity  discontinuity  and  the  face-normal  velocity  is  zero. 

The  considerable  RK2  position  error  for  the  second  of  the  two  particles  is  introduced  by  the  following  fact:  even  though 
under  perfect  time  advancement  the  second  particle  is  in  the  region  x  >  0  only  for  a  time  e  <c  At/2,  the  corrected  velocity  of 
the  RK2  step  (U(X2t0,  0)  +  U(X^\  At))/2,  is  the  average  velocity  that  the  particle  would  experience  if  it  were  in  that  region  for 
a  time  At/2.  In  other  words,  RK2  does  not  account  for  how  much  time  a  particle  spends  on  each  side  of  a  discontinuity.  This  is 
the  motivation  for  the  Multi-Step  RK2  scheme,  which  we  now  describe  in  detail. 
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2. 1.  Description  of  MRK2 

In  this  section,  we  present  a  description  of  a  Multi-Step  RK2  time  step.  We  consider  a  particle  that  is  initially  at  the  posi¬ 
tion  Xo  at  time  0,  and  which  is  to  be  advected  for  a  time  step  of  length  At.  For  the  moment,  let  us  consider  the  case  in  which 
there  is  only  one  velocity  discontinuity,  as  in  the  example  above. 

First,  we  take  an  RK2  predictor  step,  according  to  Eq.  (9): 

Xt,|  =  Xo  +  AtU(Xo,0). 

If  Xo  and  X(,)  are  in  the  same  cell,  then  we  take  the  remainder  of  a  standard  RK2  time  step,  according  to  Eq.  (10),  which  gives 
the  final  result; 

XMRK2  =  XRK2  =  ^  f  At(U(X(,  0)  +  U(X"  \  At))/2. 

Otherwise,  let  X  be  the  point  where  the  ray  from  Xo  to  X(1  first  intersects  a  cell  face.  For  Cartesian,  cylindrical  and  unstruc- 
tured  tetrahedronal  grids,  in  order  to  determine  X'  we  have  to  identify  the  faces  of  the  current  grid  cell,  each  of  which  is  a 
piece  of  a  plane  or  a  cylinder  in  2D  or  3D,  and  identify  the  points  where  the  ray  from  Xo  to  X(1)  crosses  each  of  them,  if  it  does 
(a  simple  problem  in  analytic  geometry).  Out  of  these  points,  X'  is  the  one  closest  to  Xo. 

We  make  the  following  second-order  accurate  estimate  for  the  time,  Ati,  that  it  takes  for  the  particle  to  reach  this 
discontinuity: 

At,  =At||X'-Xo||/||X")-X0||,  (11) 

where  ||  ■  ||  denotes  the  2-norm.  Next,  we  break  up  the  RK2  predictor  step  into  two  substeps.  The  first  predictor  substep,  from 
time  0  to  Ati,  and  from  Xo  to  X',  accounts  for  the  advection  of  the  particle  before  it  reaches  the  discontinuity.  The  second 
predictor  substep,  which  accounts  for  advection  after  the  particle  crosses  the  discontinuity,  takes  the  particle  from  time 
Ati  to  At  and  from  X'  to  X2),  where  X(2)  is  given  by 

X(2)  =  X'  +  (At  —  At,)lT(X'  0).  (12) 

There  are  two  possible  definitions  of  velocity  at  X',  where  the  velocity  field  is  discontinuous.  The  value  LT(X'  0)  denotes  the 
velocity  experienced  by  the  particle  after  it  has  crossed  the  discontinuity.  Similarly,  below  we  use  U  (X'  0)  to  denote  the 
velocity  experienced  by  the  particle  before  it  crosses  the  discontinuity.  (We  have  assumed,  for  the  moment,  that  the  particle 
does  not  cross  another  discontinuity  between  X'  and  X(2).) 

Next,  we  make  second-order  linear  approximations  in  time  to  estimate  the  velocities  at  X'  at  the  intermediate  time  Ati 

U  '(X\At,)«0“(X#  Ati)  =  [(At-Ati)/At]U  (X\0)  +  [Ati/At]U  (X , At),  (13) 

IT (X\  At, )  «  U+ (X',  Ati )  =  [(At  -  Ati  )/At]LT  (X'  0)  +  (At, /At]U 1  (X',  At).  (14) 

Finally,  we  calculate  the  corrected  velocities  for  each  of  the  two  substeps,  and  obtain  the  final  particle  position 

x£RK2  =Xo  +  (Af,/2)(U(Xo  0)  +  U  (X'.  At, ))  +  ((At  -  At, )/2)(0 *  (X',  At,)  +  U(X(2),  At)).  (15) 

As  already  mentioned,  for  a  particle  which  crosses  the  discontinuity,  we  break  up  the  time  step  into  two  substeps.  The 
two  terms  added  to  Xo  in  Eq.  (15)  are  respectively  the  contributions  of  each  of  those  two  substeps  to  the  advancement  of 

the  particle.  For  simplicity’s  sake,  the  algorithm  presented  above  allows  for  two  substeps  at  most,  and  can  therefore  deal 

with  at  most  one  velocity  discontinuity  crossed  by  a  particle  during  a  time  step.  By  allowing  for  a  greater  number  of  substeps 
(e.g.,  breaking  up  the  substep  from  At,  to  At  into  two  substeps  if  the  ray  from  X'  to  X(2  crosses  another  discontinuity),  the 
algorithm  is  extended  to  deal  with  an  arbitrary  number  of  discontinuities  crossed  by  a  particle  in  one  time  step. 

Below,  we  make  a  few  important  observations  about  MRK2: 

1 .  The  MRK2  procedure  is  almost  the  same  as  taking  two  RK2  steps  -  one  from  0  to  Ati  and  one  from  At,  to  At  -  with  the  one 
difference  being  that  the  MRK2  procedure  does  not  require  knowledge  of  the  velocity  field  at  the  intermediate  time  At,, 
and  uses  instead  a  linear  interpolation  in  time  between  the  velocity  fields  at  time  0  and  At.  This  difference  is  essential  in 
an  LES/FDF  context,  because  during  each  time  step  a  fraction,  approximately  equal  to  the  Courant  number,  of  the  total 
number  of  particles  in  the  domain  cross  a  cell  face.  For  a  typical  LES/FDF  calculation  with  107  particles  and  a  Courant 
number  of  0.1,  this  means  that  for  each  time  step  there  would  be  approximately  106  particles  crossing  a  cell  face,  each 
at  a  different  intermediate  time  Ati.  It  would  be  extremely  inefficient,  if  not  infeasible,  to  calculate  the  whole  recon¬ 
structed  velocity  field  at  each  of  the  (of  order  106)  intermediate  times.  Instead,  MRK2  uses  just  two  fields  to  obtain 
the  needed  local  information  with  the  required  accuracy. 

2.  For  a  fixed  grid  size,  and  subject  to  a  condition  which  is  described  in  the  next  item,  it  can  be  shown  analytically  that  an 
MRK2  step  which  crosses  a  discontinuity  introduces  an  0(At2)  error  in  the  final  particle  position,  and  preserves  the  con¬ 
tinuity  of  the  Lagrangian  mapping  X(Y,  At).  Therefore,  MRK2  preserves  continuity  (in  a  restricted  sense)  and  is  second- 
order  accurate  in  time,  since  the  number  of  discontinuities  crossed  by  a  particle  does  not  increase  with  decreasing  time 
step. 
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3.  The  condition  necessary  for  Item  2  is  that  the  velocity  normal  to  a  face  (discontinuity)  is  not  identically  zero.  MRK2  does 
not  preserve  the  continuity  of  X(Y,  At)  where  the  face-normal  velocity  is  zero.  To  see  this,  consider  a  velocity  field  similar 
to  the  example  at  the  beginning  of  Section  2,  this  time  with  D(x,y)=-0,  and  consider  two  particles:  one  at 
(Xi.o  Yi.o)  =  (—e,0)  and  the  other  at  (X2.o,  Y20)  =  (6,0).  The  result  of  an  MRK2  time  step  is  the  same  as  that  of  an  RK2  time 
step,  and  of  perfect  time  advancement:  (X1>Af,  Y i.A()  =  (-6, 0)  and  (X2iAr,  Y2&)  =  (6,  At).  Taking  limits  as  e  -*  0,  we  can  see 
why  continuity  is  violated.  Note  that  this  is  a  property  of  the  particular  velocity  field  considered,  not  a  shortcoming  of 
MRK2  (i.e.,  even  with  perfect  time  advancement  there  is  a  continuity  violation). 

4.  The  MRK2  procedure  is  applicable  to  a  wide  variety  of  grids  -  it  has  only  two  grid-dependent  aspects,  namely  the  eval¬ 
uation  of  the  interpolated  velocity  at  an  arbitrary  point  in  space,  and  the  determination  of  the  point  X'  where  the  ray  from 
Xo  to  X{1)  first  intersects  a  cell  face.  The  first  of  these  aspects,  the  ability  to  evaluate  interpolated  velocity  at  an  arbitrary 
point,  is  a  general  requirement  for  any  reasonable  velocity  interpolation  scheme.  The  second  aspect,  the  ability  to  deter¬ 
mine  efficiently  X',  depends  on  the  grid  geometry  and,  as  mentioned  in  the  description  of  M  RK2,  reduces  to  a  set  of  simple 
analytic  geometry  problems,  such  as  finding  the  intersection  of  a  ray  with  a  plane  (or  cylinder).  Therefore,  in  addition  to 
Cartesian  and  cylindrical  grids,  MRK2  in  itself  is  also  applicable  all  grids  with  planar  or  cylindrical  cell  faces  (e.g.,  unstruc¬ 
tured  tetrahedronal  grids).  However,  we  should  also  keep  in  mind  that  the  PERM  velocity  interpolation  scheme,  which  is 
the  primary  motivation  for  the  modification  from  RK2  to  MRK2,  has  currently  been  developed  only  for  Cartesian  and 
cylindrical  grids. 


3.  Comparison  between  the  performance  of  RK2  and  MRK2  when  applied  to  a  PERM  reconstructed  velocity  field 

In  this  section,  we  apply  both  RK2  and  MRK2  to  a  2D  test  problem,  and  compare  the  performance  of  the  two  time-step¬ 
ping  schemes  based  on  three  criteria:  preservation  of  continuity,  accuracy  relative  to  perfect  time  advancement,  and  the 
maintaining  of  a  consistent  subgrid  particle  distribution. 

We  use  the  following  problem  framework:  we  have  a  2D  domain,  [0, 2n\  x  [0, 2 n]t  with  periodic  boundary  conditions.  We 
consider  two  2D  periodic  flows.  Flow  1  is  incompressible,  given  by  the  formula 


U(x,y,  t)  =  1-2  cos(x  -  t)  sin(y  -  t), 

(16) 

V(x,y,  t)  =  1  +2  sin(x  -  t)  cos(y  -  t). 

(17) 

Flow  2  is  compressible,  given  by 

U(x,y.  t)  =1  -  2cos(x  -  t)sin(y  - 1)  sin(x)cos(2t). 

(18) 

V(x,y,  t)  =1  4  2  sin(x  -  t)cos(y  -  t)  4^  sin (y) cos(2t). 

(19) 

Streamline  plots  of  these  two  test  flows,  for  r  0,  are  shown  in  Fig.  2. 

In  order  to  compute  the  reconstructed  field,  the  PERM  reconstruction  scheme  uses  values  of  U  (according  to  Eq.  (16)  or 
(18))  at  the  centers  of  the  vertical  cell  faces,  and  values  of  V(according  to  Eq.  (17)  or  (19))  at  the  centers  of  the  horizontal  cell 
faces  -  this  simulates  FV  velocity  information.  Once  the  subgrid  velocity  is  reconstructed  by  PERM,  the  respective  time 
advancement  scheme  is  used  to  advect  the  particles.  The  time  step  At  is  directly  proportional  to  the  Courant  number  (C), 
which  is  defined  here  as 

C  =  max(m ax(U(x,y,  t)At/Ax)  max(V(x  y,  t)At/Ax)),  (20) 

where  Ax  is  the  side  of  a  grid  cell,  and  the  latter  two  maxima  are  over  x,  y,  t. 

3.1.  Preservation  of  continuity,  and  a  qualitative  analysis  of  particle  distributions 

Here,  we  present  results  for  the  incompressible  Flow  1,  Eqs.  (16)  and  (17).  One  of  the  beneficial  properties  of  PERM  is  that 
if  the  velocity  field  is  divergence  free  at  the  FV  level,  then  the  subgrid  reconstructed  velocity  is  also  divergence  free.  There¬ 
fore,  applying  PERM  to  Flow  1,  we  obtain  a  subgrid  velocity  field  that  is  divergence  free,  and  so  (with  perfect  time  advance¬ 
ment)  a  particle  distribution  which  is  initially  uniform  remains  uniform. 

Here,  we  advect  409,600  particles  which  are  initialized  in  the  following  manner:  the  domain  is  broken  up  into 
640  x  640  squares  of  equal  size,  and  a  single  particle  is  placed  randomly,  with  uniform  probability  density,  in  each  square. 
This  allows  for  an  initial  particle  distribution  which  is  both  very  uniform  and  does  not  degenerate  when  a  large  value  of 
strain  is  applied  to  it.  Note  that  the  number  of  particles  that  we  are  using  is  extremely  large  -  25,600  particles  per  cell,  for 
a  4  x  4  grid.  Naturally,  in  practical  LES/FDF  applications  this  number  is  much  smaller.  The  reason  why  we  are  using  so 
many  here  is  for  the  reader  to  be  able  to  better  visualize  the  mappings  X(Y,  t)  between  initial  and  final  particle  positions. 
In  the  following  figures,  each  of  the  particles  is  shown  as  a  single  light  grey  dot,  except  for  those  which  are  initially  in  the 
region  [0, 71/2]  x  [0  71/2],  at  the  lower-left  corner.  These  particles  are  dark  grey,  in  order  to  help  visualize  the  strain  of  the 
particle  distribution. 


879 8 


P.P.  Popov  et  ol./Joumol  of  Computotionol  Physics  227  (2008)  8792-8806 


Fig.  2.  Top:  streamlines,  at  t  =  0.  of  the  incompressible  Flow  1  (Eqs.  (16)  and  (17)).  Bottom:  streamlines,  at  r  =  0  of  the  compressible  Flow  2  (Eqs.  (18)  and 
(19)). 

Fig.  3  shows  results  for  a  4  x  4  grid,  with  C  1.  We  note  that  this  is  a  rather  large  Courant  number  -  it  is  used  mainly  to 
emphasize  the  effects  of  both  time  advancement  schemes.  After  a  single  time  step,  we  can  see  that  the  RK2  distribution  (top 
left  of  Fig.  3)  has  considerable  voids  where  there  are  no  particles  at  all  -  this  indicates  that  the  mapping  X(Y,  At)  is  not  sur¬ 
jective  (onto),  since  there  are  values  of  X  which  cannot  be  attained  for  any  value  of  Y.  We  can  also  see  regions  in  the  RK2 
distributions  where  the  particles  are  twice  as  dense  -  this  indicates  that  the  mapping  X(Y,  At)  is  not  injective  (one-to- 
one),  since  regions  which  were  initially  separate  are  mapped  on  top  of  each  other.  This  is  a  serious  problem  -  it  introduces 
error  in  the  subgrid  particle  distribution  (as  we  will  see  later  on)  and  violates  the  fundamental  principle  that  for  a  periodic 
flow  X(Y,  At)  is  both  injective  and  surjective. 

A  single  MRK2  time  step,  too,  has  these  problems,  but  to  a  much  lesser  extent  -  we  can  see  (from  the  top  right  of  Fig.  3) 
four  elongated  gaps  where  there  are  no  particles,  but  the  size  of  these  gaps  is  much  smaller  than  the  size  of  the  holes  in  the 
particle  distribution  created  by  RK2.  These  gaps  in  the  MRK2  distribution  are  to  be  expected  -  they  correspond  to  regions 
where  the  face-normal  velocity  is  zero,  as  explained  in  Section  2.  Comparing  the  particle  distributions  after  one  flow-through 
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Fig.  3.  Particle  distributions  in  Flow  1  for  a  4  x  4  grid,  C  -  1 .  Top  left:  RK2.  one  time  step  (At  =  0.448);  top  right:  MRK2.  one  time  step;  bottom  left.  RK2,  one 
flow-through  time  (12  time  steps);  bottom  right:  MRK2.  one  flow-through  time.  The  dark  grey  particles  are  those  that  are  in  the  lower-left  corner  of  the 
domain  at  t  -  0. 


time  (12  time  steps),  we  can  see  the  same  differences  but  with  greater  magnitude.  Whereas,  MRK2  does  not  preserve  the 
uniformity  of  the  particle  distribution  very  well,  it  does  much  better  than  RK2,  where  there  are  considerably  larger  void 
regions. 

Fig.  4  shows  results  for  the  same  value  of  the  Courant  number,  C  —  1 ,  but  for  a  finer  8x8  grid.  Because  the  magnitudes  of 
the  discontinuities  in  PERM  are  0(Ax2)  [1  ],  the  velocity  discontinuities  in  this  case  are  theoretically  smaller  by  a  factor  of 
four,  and  RK2  does  almost  as  well  as  MRK2. 

Fig.  5  shows  results  for  a  4  x  4  grid,  with  a  smaller  Courant  number,  C  0  25.  We  can  again  see  that  RK2  does  not  pre¬ 
serve  continuity  very  well  -  there  are  noticeable  voids  in  the  RK2  distribution,  both  after  one  time  step  and  after  one  flow¬ 
through  time.  Although  we  know  theoretically  that  there  are  also  voids  in  the  MRK2  distributions,  these  cannot  be  perceived 
on  the  plots.  We  also  note  that  MRK2  performs  better  at  preserving  the  sharp  interface  between  the  light  and  dark  grey 
regions. 

Fig.  6  shows  results  for  an  8  x  8  grid,  C  =  0  25.  We  cannot  perceive  a  difference  between  the  two  particle  distributions, 
which  suggests  that,  due  to  the  smaller  value  of  velocity  discontinuities,  the  performance  of  RK2  is  similar  to  that  of  MRK2. 

From  these  comparisons  of  particle  distributions  produced  by  RK2  and  MRK2,  we  conclude  that  when  the  velocity  discon¬ 
tinuities  yielded  by  the  reconstruction  are  considerable  (as  for  the  4x4  grid),  MRK2  performs  better  than  RK2  at  preserving 
the  continuity  of  the  particle  distribution  and  the  consistency  between  the  subgrid  particle  distribution  and  the  LES  filtered 
density. 

3.2.  Accuracy  relative  to  perfect  time  advancement 

As  we  are  focused  on  developing  and  demonstrating  an  effective  time  advancement  scheme  for  a  given  velocity  recon¬ 
struction,  we  are  interested  in  the  error  of  a  given  solution  relative  to  perfect  time  advancement.  Since  we  do  not  have 
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Fig.  4.  Particle  distributions  in  Flow  1  for  an  8  x  8  grid.  C  « 1.  Top  left:  RK2,  one  time  step  (At  -  0.224),  top  right  MRK2,  one  time  step;  bottom  left:  RK2. 
one  flow-through  time  (24  time  steps);  bottom  right;  MRK2.  one  flow-through  time.  The  dark  grey  particles  are  those  that  are  in  the  lower-left  corner  of  the 
domain  at  f  =  0. 


an  exact  perfect  time  advancement  solution,  we  use  instead  an  MRK2  solution  with  a  very  low  Courant  number  (C  1/512) 
as  an  approximation. 

First  of  all,  we  consider  position  error.  We  use  N  256  particles,  uniformly  distributed  across  the  domain  at  time  0.  We 
denote  by  Y(0  the  initial  position  of  the  ith  particle,  and  we  denote  by  XRK2(Y0),  t),XMRK2(Y(l),  t),  and  X(Y(,),r)  the  correspond¬ 
ing  Lagrangian  position  mappings  yielded  by  RK2,  MRK2  and  perfect  time  advancement,  respectively.  The  position  error  of  a 
given  RK2  solution  is  then  defined  as 

e*RK2  =  |Xriq(Y(,),  t)  -  X(Y"V)ll)-  (21) 

The  MRK2  position  error  is  defined  analogously.  As  we  already  mentioned,  it  can  be  shown  analytically  that,  with  At  being 
the  length  of  the  time  step,  the  RK2  position  error  is  of  order  O(At),  and  the  MRK2  position  error  is  of  order  0(At2). 

We  also  consider  error  of  infinitesimal  volume  expansion  (which  we  will  refer  to  as  d  Verror,  for  the  sake  of  brevity).  Note 
that  we  refer  here  to  two-dimensional  volume,  i.e.,  area.  For  an  infinitesimal  material  element  whose  initial  volume  is  dV'o, 
and  whose  initial  position  is  Y,  the  final  volume,  dVtf  is  given  by 

dV,  =  dV0  det(^™) ,  (22) 

where  X(Y,  t)  is  the  Lagrangian  position  mapping  previously  mentioned.  This  is  the  basis  of  the  change  of  variables  formula, 
which  in  this  context  tells  us  that  the  volume  at  time  t  of  a  set  S  is 


(23) 
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Fig.  5.  Particle  distributions  in  Flow  1  for  a  4  x  4  grid,  C  -  0.25.  Top  left  RK2,  one  time  step  (At  -  0  112);  top  right  MRK2.  one  time  step;  bottom  left:  RK2, 
one  flow-through  time  (48  time  steps);  bottom  right:  MRK2,  one  flow-through  time.  The  dark  grey  particles  are  those  that  are  in  the  lower-left  corner  of  the 
domain  at  t  =  0. 


provided  that  the  mapping  X(Y  t)  is  one-to-one.  Therefore,  accurate  values  of  det(flX;/0Yk)  are  essential  for  achieving  an 
accurate  subgnd  particle  distribution.  We  also  note  that  X(Y,  t)  may  not  be  one-to-one  if  its  continuity  is  violated  (such 
as  in  Fig.  3,  where  particles  are  twice  as  dense  in  some  places,  because  certain  regions  are  mapped  on  top  of  each  other), 
and  therefore  continuity  of  the  position  mapping,  which  is  considered  in  the  previous  subsection,  is  just  as  essential  for  get¬ 
ting  the  right  value  of  Vt  as  are  accurate  values  of  det(0X//0Yk). 

With  this  in  mind,  we  define  dV  error  of  a  given  RK2  solution  as 


e 


RK2 

dV 


1  N 


det 


-  det 


(24) 


The  MRK2  dV  error  is  similarly  defined.  For  RK2,  it  can  be  proven  analytically  that  the  dV  error  is  0(At2)  when  the  divergence 
of  the  PERM  reconstructed  field  is  continuous  everywhere,  and  O(At)  when  the  divergence  is  discontinuous  from  cell  to  cell, 
which  is  the  case  for  all  compressible  flows. 

Here,  we  cannot  calculate  the  necessary  Jacobians  exactly,  but  a  sufficiently  accurate  numerical  estimate  which  we  use  is 

det®)«V*0,  (25) 

where  Ao  is  the  initial  area  of  a  triangle  initially  of  sides  10  8  formed  by  three  particles  at  time  0,  and  At  is  the  area  of  the 
triangle  formed  by  those  three  particles  at  time  t. 

Here,  we  present  results  for  the  compressible  Flow  2,  Eqs.  (18)  and  (19).  Fig.  7  shows  plots  of  position  error  and  dV  error 
for  4  x  4  and  8x8  grid  solutions,  with  the  Courant  number  ranging  from  1  down  to  1/32.  For  both  position  error  and  dV 
error,  RK2  exhibits  first-order  behavior,  whereas  MRK2  exhibits  second-order  behavior.  This  is  consistent  with  previously 
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Fig.  6.  Particle  distributions  in  Flow  1  for  an  8  x  8  grid.  C  -  0.25.  Top  left.  RK2,  one  time  step  ( At  =  0.056);  top  right;  MRK2,  one  time  step;  bottom  left  RK2 
one  flow-through  time  (96  time  steps);  bottom  right;  MRK2.  one  flow-through  time  The  dark  grey  particles  are  those  that  are  in  the  lower-left  corner  of  the 
domain  at  t  =  0. 


stated  analytic  results,  and  it  also  indicates  that  for  MRK2,  dV error  is  0(At2)  -  a  result  which  we  have  not  been  able  to  prove 
rigorously. 

It  should  also  be  noted  that,  for  larger  values  of  C,  MRK2  is  computationally  slower  than  RK2.  In  order  to  compare  the 
efficiency  of  MRK2  with  that  of  RK2,  Fig.  8  shows  plots  of  the  position  and  d^  errors  from  Fig.  7  versus  simulation  time, 
for  a  simulation  with  1600  particles  per  cell,  instead  of  just  the  256  particles  (for  the  entire  domain)  that  we  used  to  deter¬ 
mine  the  errors.  We  use  a  larger  number  of  particles  here  to  ensure  that  the  cost  of  all  the  operations  other  than  particle 
advection  (such  as  the  determination  of  the  PERM  coefficients)  is  negligible. 

From  Fig.  8  we  see  that  MRK2  becomes  more  efficient  as  the  error  (and  hence  the  Courant  number)  decreases  -  for  the 
4x4  grid,  MRK2  is  more  efficient  for  C  O  /2,  and  for  the  8  x  8  grid,  MRK2  is  more  efficient  for  C  ^  1  /4.  On  the  other  hand, 
Fig.  8  takes  into  account  only  the  computational  cost  for  particle  advection  -  if  we  also  consider  the  cost  of  chemistry  cal¬ 
culations,  which  is  usually  much  higher  than  the  cost  of  particle  advection,  we  conclude  that  MRK2  may  be  more  efficient 
than  RK2  even  for  larger  values  of  C. 

Table  1  shows  values  of  the  error,  as  well  as  simulation  time  for  a  1600-particles-per-cell  simulation,  for  C  1  and  for 
C  -  0.25.  For  both  values  of  C,  and  for  both  meshes,  MRK2  has  a  smaller  position  error:  f^RK2/eJK2  <  0.69.  For  both  meshes, 
the  MRK2  dV  error  is  less  than  the  RK2  dV error  for  C  =  0.25,  but  for  C  =  1  and  the  8  x  8  mesh,  the  MRK2  d^  error  is  greater 
than  the  RK2  dV  error.  This,  however,  is  not  a  major  concern,  because  as  we  already  noted,  C  =  1  is  a  very  large  value  for  the 
Courant  number  -  in  most  applications,  we  expect  a  smaller  value,  such  as  0.25,  to  be  used. 

3.3.  Accuracy  of  the  subgnd  particle  distribution 

Previously,  McDermott  and  Pope  [  1  ]  have  shown  that  a  PERM  reconstructed  velocity  field,  with  RK2  as  the  time-stepping 
scheme,  performs  much  better  than  linear  interpolation  in  maintaining  a  subgnd  particle  distribution  which  is  consistent 
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Fig.  7.  Results  for  position  and  d terror  in  Flow  2  at  r  =  8.97.  Top  left,  position  error  for  a  4  x  4  grid;  top  right:  dV error  for  a  4  x  4  grid;  bottom  left,  position 
error  for  a  8  x  8  grid;  bottom  right;  dV  error  for  a  8  x  8  grid  The  diamonds  denote  RK2  error,  the  circles  denote  MRK2  error.  The  dashed  line  indicates  first- 
order  convergence,  and  the  solid  line  indicates  second-order  convergence. 


with  the  LES  filtered  density.  Here,  we  demonstrate  that  this  consistency  becomes  even  better  when  MRK2  is  used  as  the 
time-stepping  scheme. 

For  the  present  tests,  we  use  the  incompressible  Flow  1,  Eqs.  (16)  and  (17).  The  particles  are  initialized  in  a  manner  similar 
to  that  in  Section  3.1  (by  breaking  up  the  domain  into  as  many  squares  as  there  are  particles,  and  initializing  a  single  particle 
randomly  into  each  square).  Then  the  mean  particle  mass  density  is  initially  the  same  everywhere  in  the  flow,  and  it  remains 
so  under  perfect  time  advancement,  for  a  divergence-free  velocity  field  such  as  the  one  we  consider.  Hence,  the  expected 
number  of  particles  in  a  set  S  is  directly  proportional  to  the  area  of  S. 

We  quantify  the  accuracy  of  the  subgrid  particle  distribution  in  the  following  way:  we  break  up  the  domain  into  equal 
sized  squares,  which  we  call  sampling  squares.  We  use  sampling  squares  of  three  different  sizes,  with  sides  71/4,  7t/8,  and 
tt/1  6  (respectively  1  /2, 1  /4  and  1  / 8  of  the  sidelength  of  a  grid  cell,  for  a  4  x  4  grid).  At  any  given  timestep,  let  N,  be  the  num¬ 
ber  of  particles  in  the  ith  sampling  square,  let  Nm„n  be  the  expected  number  of  particles  for  a  square  of  this  size,  and  let  Nsq 
be  the  number  of  sampling  squares  in  the  domain.  The  subgrid  particle  distribution  error  is  then  defined  as 

|  ( Ni  —  Nmcan)/Wm<ran  l^/^sq-  (26) 

This  subgrid  distribution  error  has  two  contributions:  first,  a  probabilistic  error,  due  to  the  random  particle  initialization: 
and,  second,  a  bias  error,  due  to  the  particle  advection  scheme  (an  example  of  this  are  the  considerable  voids  in  the  particle 
distribution  left  by  RK2  at  C  =  1).  For  this  test  flow,  the  PERM  reconstructed  velocity  field  has  zero  divergence,  and  therefore 
a  perfect  time  advancement  solution  contains  only  the  probabilistic  component  of  the  error.  Therefore,  we  compare  the  sub- 
grid  distribution  errors  given  by  the  MRK2,  RK2  and  perfect  time  advancement  solutions,  making  the  assumption  that  what¬ 
ever  is  in  excess  of  the  perfect  time  advancement  error  is  mostly  bias  error. 

Fig.  9  shows  results  for  65,536  particles,  C  =  0.25  and  a  4  x  4  grid.  As  may  be  seen,  the  MRK2  subgrid  distribution  error  is 
much  closer  to  the  perfect  time  advancement  error  than  is  the  RK2  error.  From  this  we  deduce  that,  for  the  coarse  4x4  grid, 
the  bias  error  caused  by  MRK2  is  smaller  than  the  bias  error  caused  by  RK2. 
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Fig.  8.  Comparison  between  the  efficiency  of  MRK2  and  RK2,  based  on  position  and  dV  error,  and  the  computational  cost  of  particle  advection  for  a  1600- 
particles-per-cell  simulation  of  Flow  2  at  r  =  8.97.  Top  left:  position  error  for  a  4  x  4  grid;  top  right:  dVerror  for  a  4  x  4  grid;  bottom  left:  position  error  for  a 
8x8  grid;  bottom  right:  dV  error  for  a  8  x  8  grid.  The  diamonds  denote  RK2  error,  the  circles  denote  MRK2  error.  The  dashed  line  indicates  first-order 
convergence,  and  the  solid  line  indicates  second-order  convergence. 


Table  1 

Summary  of  position  and  dV  errors  for  Flow  2  at  f  =  8  97 


Case 

€x 

*dv 

RK2 

MRK2 

RK2 

MRK2 

4  x  4,  C  =  1 

1.32  x  10° 

8.17  x  10"1 

4.35  x  10  1 

3.97  x  10"1 

4  x  4.  C  =  0  25 

4.35  x  10  1 

1 .29  x  10" 1 

1.07  x  10  1 

4.56  x  10~2 

8  x  8,  C  =  1 

5.47  x  10  1 

3.80  x  10  1 

1.88x10  1 

2  33  x  10  1 

8  x  8,  C  =  0  25 

7.83  x  10'2 

3.55  x  10"2 

1.84x10  2 

1  41  x  10  2 

Fig.  1 0  shows  results  for  65,536  particles,  C  0.5  and  a  8  x  8  grid.  The  time  step  is  the  same  as  for  the  results  on  Fig.  6,  but 
the  velocity  field  is  better  resolved  by  PERM.  Here,  the  MRK2,  RK2  and  perfect  time  advancement  errors  are  much  closer  to¬ 
gether,  and  MRK2  does  not  provide  a  big  advantage  over  RK2.  Again,  we  observe  that  for  a  finer  grid  RK2  does  not  perform 
considerably  worse  than  MRK2,  because  the  discontinuities  in  velocity  are  smaller. 

3.4.  Computational  cost 

We  have  seen  that  MRK2  has  superior  accuracy  than  RK2,  and  has  better  performance  in  preserving  the  consistency  be¬ 
tween  the  subgrid  particle  distribution  and  the  filtered  LES  density.  One  disadvantage  that  MRK2  has  over  RK2  is  that  it  is 
slower  -  this  was  mentioned  in  Section  3.2,  where  we  compared  the  efficiency  of  MRK2  with  that  of  RK2. 

For  the  well-optimized  Fortran  90  implementation  of  MRK2  which  we  used  to  obtain  the  results  in  the  previous  sections, 
each  additional  substep  after  the  first  one  costs  approximately  1.5  times  the  cost  of  an  RK2  step.  In  other  words,  if  the  total 
time  for  one  RK2  time  step  for  one  particle  is  tstep,  the  total  time  for  an  MRK2  time  step,  for  a  particle  which  crosses  n  dis¬ 
continuities,  is  approximately  rstep(l  +  1.5n).  Therefore,  for  a  flow  in  which  each  particle  crosses  one  discontinuity  per  time 
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Fig.  9.  Results  for  subgrid  particle  distribution  error  for  Flow  1.  C  -  0.25,  (Ar  =  0.1 12).  4x4  grid  top:  error  for  sampling  squares  of  side  ti/4;  middle:  error 
for  sampling  squares  of  side  tt/8;  bottom:  error  for  sampling  squares  of  side  ti/16.  The  errors  for  the  RK2,  perfect  time  advancement,  and  MRK2  solutions 
are  shown  respectively  by  the  dashed,  solid,  and  dash-dot  curves. 


Fig.  10.  Results  for  subgrid  particle  distribution  error  for  Flow  1.  O  0.5,  (At  =  0.1 12).  8x8  grid  top:  error  for  sampling  squares  of  side  ti/4;  middle*  error 
for  sampling  squares  of  side  7T/8;  bottom:  error  for  sampling  squares  of  side  ti/16.  The  errors  for  the  RK2,  perfect  time  advancement,  and  MRK2  solutions 
are  shown  respectively  by  the  dashed,  solid,  and  dash-dot  curves. 


step  (hence  requiring  two  substeps  per  time  step)  an  MRK2  calculation  would  require  2.5  times  more  time  than  an  RK2  cal¬ 
culation  would. 

Fortunately,  in  most  flows  of  practical  interest  the  number  of  particles  which  require  two  or  more  substeps  is  generally 
small.  For  the  test  cases  from  the  previous  sections,  the  MRK2  calculations  have  a  70%  computational  overhead  at  C  1 ,  and 
a  20%  computational  overhead  at  C  =  0.25,  relative  to  the  RK2  calculations.  For  smaller  values  of  C  the  computational  over¬ 
head  is  even  smaller,  since  the  ratio  between  the  number  of  particles  which  cross  a  discontinuity  and  the  number  of  particles 
which  do  not  cross  one  decreases  with  decreasing  Courant  number. 

4.  Conclusions 

We  have  outlined  a  new  time-stepping  scheme,  based  on  second-order  Runge-Kutta,  which  is  particularly  suited  for 
advection  of  a  large  number  of  particles  through  a  discontinuous  velocity  field,  such  as  the  one  yielded  by  the  parabolic  edge 
reconstruction  method  (PERM).  We  have  demonstrated  that  the  new  scheme,  MRK2,  preserves  better  than  RK2  the  continu¬ 
ity  of  the  Lagrangian  mapping  between  initial  and  final  positions.  We  have  also  shown  that  MRK2  is  second-order  accurate  in 
time,  as  opposed  to  RK2,  which  is  first-order  accurate  for  a  discontinuous  velocity  field.  Additionally,  we  have  shown  that 
MRK2  preserves  better  than  RK2  the  consistency  between  the  subgrid  particle  distribution  and  the  LES  filtered  density 
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On  the  other  hand,  we  have  seen  that  the  advantages  of  MRK2  over  RK2  diminish  when  a  more  refined  grid  is  used  for 
velocity  reconstruction,  and  that  MRK2  has  a  computational  overhead  relative  to  RK2  (20%  at  C  —  0  25)  which  is  not  negli¬ 
gible.  Therefore,  whereas  RK2  is  sufficiently  accurate  (and  faster)  for  problems  with  a  fine  mesh,  it  is  preferable  to  use  MRK2 
in  problems  where  the  mesh  provides  only  marginal  resolution  of  the  velocity,  which  is  often  the  case  in  LES/FDF  methods. 

Finally,  it  should  be  noted  that  the  main  principle  of  MRK2  -  breaking  up  a  step  into  two  or  more  substeps  every  time  a 
velocity  discontinuity  is  encountered  -  can  be  applied  in  a  straightforward  manner  to  time  advancement  schemes  other  than 
RK2.  The  Midpoint  Euler  scheme,  for  example,  can  be  modified  to  deal  with  velocity  discontinuities  similarly  to  the  way 
MRK2  does.  This  would  be  preferable  in  a  situation  where  velocity  information  is  known  at  the  middle  of  the  time  step. 
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